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ABSTRACT 


The  problem  of  reconstructing  an  unknown  image  from  partial  data  in 
the  Fresnel  zone  region  is  discussed.  The  Fresnel  region  of  a  coherent 
optical  processing  system  is  used  as  a  basis  for  the  definition  of  a  Fresnel 
zone  transform  pair.  An  iterative  method  of  reconstructing  an  unknown 
image  located  in  the  object  domain,  from  either  the  phase-only  or 
magnitude-only  data  in  the  transform  domain,  is  investigated  through  the 
use  of  digital  calculations. 

The  Fresnel  zone  transform  is  implemented  digitally  using  a  fast 
Fourier  transform  routine.  The  number  of  Fresnel  zones,  which  is  a 
measure  of  how  far  out  of  the  Fourier  transform  plane  one  is  located,  is 
limited  by  certain  sampling  considerations  which  are  discussed.  From  the 
Fourier  transform  plane  through  the  Fresnel  region,  it  is  demonstrated 
that  the  phase-only  information  in  the  transform  allows  for  a  very  good 
reconstruction  in  just  a  few  iterations  of  the  reconstruction  algorithm. 
While  reconstructing  an  image  from  the  magnitude  of  the  Fourier 
transform  is  very  difficult,  as  one  moves  further  into  the  Fresnel  zone 
region,  good  reconstructions  are  obtained  in  only  several  iterations. 

The  success  of  the  magnitude-only  reconstructions  from  data  far  into 
the  Fresnel  region  is  explained.  An  approximation  of  the  Fresnel  zone 
transform  integral  is  derived  which  is  valid  for  a  large  number  of  Fresnel 
zones.  This  approximation  is  used  to  show  that  the  first  iterative  guess  is 
approximately  equal  to  the  unknown  image. 


Phase-only  information  in  the  Fresnel  zone  transform  leads  to  very 
good  reconstructions  throughout  the  Fresnel  zone  region.  For  a  given 
number  of  iterations,  there  is  a  slight  degradation  in  the  quality  of  the 
reconstructed  image  as  the  number  of  Fresnel  zones  is  increased.  The 
phase  of  the  Fresnel  zone  transform  preserves  the  structure  of  the  sharp 
edges  in  the  image,  the  ability  of  the  phase  information  to  preserve  the 
edge  information  is  demonstrated  analytically. 

A  degradation  in  the  reconstructed  image  results  if  noise  is  present  in 
the  given  data.  The  probability  density  functions  of  the  magnitude  and 
phase  of  a  signal  in  the  presence  of  coherently  additive  noise  are  derived. 
The  effects  of  both  signal  independent  and  signal  dependent  noise  are 
investigated. 

An  error  in  the  number  of  Fresnel  zones  is  investigated.  A  method  of 
searching  through  the  Fresnel  region  to  find  the  correct  location  of  the 
Fresnel  zone  transform  is  discussed. 

In  reconstructing  an  image  from  either  the  phase-only  or  magnitude- 
only  data  of  the  Fresnel  zone  transform,  the  iterative  reconstruction 
algorithms  which  are  used  are  shown  to  be  stable  both  in  the  presence  of 
noise  and  in  an  error  in  the  number  of  Fresnel  zones.  In  addition,  while 
there  is  only  a  slight  degradation  in  the  quality  of  the  reconstruction  from 
phase-only  data,  there  is  a  dramatic  improvement  in  the  quality  of  the 
magnitude-only  reconstruction  as  one  moves  out  of  the  Fourier  transform 
plane  and  into  the  Fresnel  region. 

The  novel  aspect  of  this  research  is  the  extension  of  the  reconstruction 
algorithm  into  the  Fresnel  zone  region.  With  the  formulation  presented 


VI 


in  the  appendix,  the  Fresnel  zone  transform  is  seen  to  be  a  particular  case  of  a  more 
general  reconstruction  problem.  It  is  demonstrated  in  this  dissertation  that  for  a  set 
number  of  iterations,  the  reconstructions  are  much  improved  if  the  illuminating 
wavefront  is  known  and  non-trivial  (i.e.  contains  terms  of  higher  order  than 
constant  and  linear). 
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Chapter  1 


Introduction 

1 .1  Statement  of  the  Problem  and  Objectives 

Down  through  history,  people  have  been  trying  to  get  something  for 
nothing.  It  may  be  as  simple  as  a  free  lunch  or  as  complicated  as  the 
reconstruction  of  an  unknown  image  from  partial  information  of  the 
transform.  It  is  worthwhile  to  begin  this  discussion  with  an  explanation  of 
image  reconstruction  and  how  it  relates  to  other  fields  such  as  image 
enhancement,  image  restoration,  image  synthesis,  and  image  recovery. 
The  fact  that  all  of  these  areas  are  of  active  interest  is  affirmed  by  the 
recent  publication  of  two  books  dedicated  to  these  subjects  [Stark,  (1987); 
Bates  and  McDonnell,  (1986)],  and  the  success  of  recent  conferences  [see 
for  example.  Signal  Recovery  and  Synthesis  with  Incomplete  Information 
and  Constraints  Technical  Digest,  (Optical  Society  of  America,  Washington 
D.C.,  1983);  Topical  Meeting  on  Signal  Recovery  and  Synthesis  II  Technical 
Digest,  (Optical  Society  of  America,  Washington  D.C.,  1986);  SPSE's  40th 
Annual  Conference  and  Symposium  on  Hybrid  Imaging  Systems,  (SPSE  - 
The  Society  for  Imaging  Science  and  Technology,  Springfield  VA,  1987)] 
and  special  issue  journal  publications  dedicated  to  these  topics  [see  for 
example,  J.  Opt.  Soc.  Am.,  73,  (1983);  J.  Opt.  Soc.  Am.  A ,  4,  (1987)]. 

Manipulating  the  tone  reproduction  curve,  adjusting  color  maps  for 
the  display  of  a  digital  image,  and  changing  the  development  process 
when  printing  a  picture  are  just  some  examples  of  image  enhancement. 


"Image  enhancement  processes  consist  of  a  collection  of  techniques  that 
seek  to  improve  the  visual  appearance  of  an  image,  or  to  convert  the 
image  to  a  form  better  suited  to  human  or  machine  analysis*  [Pratt, 
(1978),  p.  307].  Many  image  enhancement  techniques  lack  a  rigorous 
derivation  as  a  basis  for  the  methods  which  are  used.  This  is  quite  often 
the  case  because  the  end  result,  the  improved  appearance  of  an  image,  is 
very  subjective. 

Image  restoration  is  very  simitar  to  image  enhancement,  in  that  the 
result  is  an  improved  image.  "Image  restoration  involves  estimating  the 
properties  of  the  distortion  and  using  them  to  refurbish  the  original 
information"  [Bates  and  McDonnell,  (1986),  p.  5].  The  techniques  used  in 
image  restoration  generally  have  a  mathematical  basis  and  involve 
estimation,  statistics,  and  the  modeling  of  physical  processes.  Some 
typical  image  restoration  procedures  involve  estimating  the  statistical 
properties  of  various  noise  processes  in  an  attempt  to  read  a  signal 
through  the  noise,  or  the  modeling  of  the  physical  effects  of  varnish  and 
shellac  over  time  to  restore  ancient  artwork  [Asmus,  (1987)]. 

In  both  the  enhancement  and  restoration  of  images,  one  is  given  an 
image  and  work  is  done  to  improve  that  image  for  its  intended  use. 
Image  synthesis  differs  from  these  procedures  because  there  is  no  original 
image  given;  "the  synthesis  problem  is  concerned  with  the  determination 
of  a  signal,  such  as  an  image,  which  is  consistent  with  the  known 
constraints"  [Hayes,  (1987)].  "An  image  restoration  problem  involves  the 
recovery  of  an  actual  image  that  has  been  distorted  by  a  given  imaging 
system.  Image  synthesis,  on  the  other  hand,  aims  at  discovering  an  input 


image  which,  when  distorted  by  a  given  imaging  system  (e.g.  a  display 
device),  yields  a  desired  output  pattern.  The  sought  image  is  discovered 
(designed,  constructed,  or  synthesized),  rather  than  recovered, 
reconstructed,  or  restored"  [Saleh,  (1987)].  A  direct  application  of  these 
synthesis  techniques  is  in  the  design  of  a  mask  used  for  the  manufacture 
of  integrated  circuits.  In  this  instance,  there  is  a  very  definite  structure 
desired  after  the  imaging  of  a  mask  onto  a  substrate,  and  the  subsequent 
etching  of  that  substrate. 

"Stated  in  its  most  general  form,  the  signal  recovery  problem  is 
described  like  this:  Given  data  g,  determine  the  source  f,  that  produced  g" 
[Stark,  (1987),  p.  xi].  This  definition  of  image  recovery  is  indeed  very 
general,  and  may  better  be  viewed  as  just  one  topic  under  the  wider  field 
of  inverse  source  problems.  In  the  image  enhancement  and  image 
restoration  problems  discussed  above,  one  is  given  an  image  of  some  type, 
and  the  procedures  will  alter  the  image  in  some  way,  but  the  general 
structure  and  form  of  the  image  is  known  and  not  changed.  In  the  case  of 
image  recovery,  the  given  data  may  differ  quite  a  bit  from  the  sought 
after  solution  to  the  problem.  Literature  concerning  both  the  larger  areas 
of  study  investigating  various  aspects  of  of  the  inverse  source  problem, 
and  inverse  scattering  problem  (see  for  example:  Baltes,  ed.,  (1978); 
Baltes,  ed.,  (1980)],  and  the  general  field  of  image  recovery  will  not  be 
reviewed  here. 

The  work  contained  in  this  thesis  is  best  placed  in  the  area  of  image 
reconstruction.  "The  reconstruction  problem  ...  is  concerned  with  the 
unique  recovery  of  a  signal  given  a  partial  Fourier  domain  description 
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along  with  a  signal  model,  or  a  set  of  signal  constraints*  {Hayes,  (1987)]. 
The  reconstruction  of  a  signal  from  partial  Fourier  domain  information 
includes  the  areas  of  computer  aided  tomography,  image  retrieval,  phase 
retrieval,  and  magnitude  retrieval.  The  use  of  Fresnel  zone  information  as 
opposed  to  Fourier  domain  information  is  a  novel  aspect  of  the  research 
discussed  in  this  dissertation. 

While  the  Fresnel  zone  transform  pair  will  be  defined  in  Sect.  2.2,  it  is 
worthwhile  to  define  the  Fourier  transform  pair.  Given  a  function  (or 
image),  g(x,y),  defined  in  the  object  domain  (or  image  plane),  the  Fourier 
transform,  G(fx,  fy),  is  given  by 

-•2n(/,x  +  r.y)  .  .  (1-1) 


G(^^,/y)=  j  j  g(x^)e  ‘  ^  dxdy  , 

where  x  and  y  are  spatial  coordinates  and  fx  and  fy  are  spatial  frequencies. 
The  inverse  Fourier  transform  is  given  by 

+i2n{fx  +  fy)  (1-2) 

«(*.y)=  J  e  '  ^  dfdr^.  '  ' 

Throughout  this  thesis,  lower  case  letters  will  be  used  to  represent 
functions  (or  images)  in  the  object  domain  (or  image  plane),  and  upper 
case  letters  will  represent  the  transform  of  the  corresponding  lower  case 
function.  Writing  the  complex  valued  function,  G(fx,fy),  in  polar  form 


=  .  (1-3) 

the  magnitude  (or  modulus),  |G(fx,  fy)|,  is  given  by 

=  ,  (1-4) 

where  *  represents  a  conjugate  complex.  The  phase,  e(fx,  fy),  is  defined  to 
satisfy  the  simultaneous  equations 


(1-5a) 


\Gif^.f^)\ 


(1-5b) 


where  RE[G(fx.fy)]  and  IM[G(fx.fy)]  are  the  real  and  imaginary  parts, 
respectively,  of  G(fx,fy). 

The  terms  magnitude-only  reconstruction  and  phase  retrieval  are  both 
used  to  describe  the  problem  when  only  the  magnitude  information, 
|G(fx,fy)|.  is  known.  The  problem  is  then  one  of  reconstructing  the  image 
.'rom  this  magnitude-only  information,  or  retrieving  the  the  phase  which 
corresponds  to  the  given  magnitude.  The  terms  phase-only  reconstruction 
and  magnitude  retrieval  are  both  used  to  describe  the  problem  when  only 
the  complex  phase,  exp{  I  e(fx,  fy) },  is  known.  The  problem  is  then  one  of 
reconstructing  the  image  from  this  phase-only  Information,  or  retrieving 
the  magnitude  which  corresponds  to  the  given  phase  information.  A  third 
type  of  problem  is  posed  when  both  the  magnitude  and  phase  are  known, 
but  only  over  a  partial  region  of  the  transform  space.  This  is  the  problem 
often  encountered  in  computer  aided  tomography.  All  of  these  problems 
generally  make  certain  assumptions  about  the  unknown  image  such  as 
constraints  on  its  spatial  extent  or  non-negativity. 


1.2  Background 

The  problem  of  reconstructing  a  function  from  partial  Fourier  domain 
information  was  first  investigated  by  A.A.  Michelson  [1891,  1892]  and 
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Lord  Rayleigh  [1892].  The  problem  originally  arose  in  trying  to  determine 
the  spectrum  of  a  source  from  a  measurement  of  the  visibility  curves  of 
the  interference  pattern.  The  measured  visibility  curve  represented  only 
the  modulus  of  the  Fourier  transform  of  the  unknown  spectrum  [see  also 
Born  and  Wolf.  §7.5.8,  (1980)]. 

The  phase  retrieval  problem  was  investigated  by  Wolf  [1962],  who 
introduced  the  idea  of  using  some  of  the  analytic  properties  of  the 
functions  involved.  This  method  used  the  locations  of  the  complex  zeros 
in  the  analytically  continued  functions  to  reconstruct  the  complete  Fourier 
spectra. 

The  dispersion  relations  [Toll,  (1956)],  which  relate  the  real  and 
imaginary  parts  of  a  complex  function  that  is  analytic  and  regular  in  the 
lower  half  plane,  are  used  to  reconstruct  the  Fourier  phase  from 
knowledge  of  the  log  of  the  Fourier  magnitude.  Problems  arise  however, 
when  there  are  zeros  located  in  the  lower  half  plane  of  the  analytically 
continued  Fourier  magnitude.  In  this  case  there  is  a  multiplicity  of 
solutions  due  to  the  Blaschke  factor  which  is  used  to  remove  these  zeros. 
The  position  of  the  zeros  and  their  physical  interpretation  was 
investigated  by  Nussenvieg  [1966].  In  addition,  it  has  been  shown  [Nieto- 
Vesperinas,  (1980)]  that  the  extension  of  the  dispersion  relations  to  two 
dimensions  yields  a  pair  of  equations  that  in  general  have  no  solution. 

Techniques  which  made  use  of  the  fact  that  the  functions  are  band- 
limited,  which  is  the  case  for  imaging,  were  investigated  by  Walther 
[1962].  This  method  dealt  with  only  the  one  dimensional  diffraction 
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problem,  and  Walther  expected  that  "the  extension  to  the  real  two- 
dimensional  case  will  probably  be  fraught  with  difficulties." 

Methods  of  modifying  the  experimental  measurement  were  first 
suggested  by  Kohler  and  Mandel  [1969,  1972],  and  Mehta  [1965,  1968]. 
These  methods  suggested  the  use  of  an  additional  point  source  or 
reference  beam,  which  is  analogous  to  the  use  of  a  reference  beam  in 
holography,  and  the  use  of  exponential  filters  on  the  source,  which  would 
change  the  locations  of  the  zeros  in  the  analytically  continued  complex 
plane.  The  conditions  under  which  the  support  of  the  unknown  object 
allow  for  the  successful  reconstruction  were  also  investigated 
[Greenaway,  (1977)]. 

Using  the  location  of  the  zeros  not  only  in  the  Fourier  transform 
magnitude,  but  in  the  Fresnel  region  as  well,  has  also  been  suggested 
[Ross  et  al  (1980)],  along  with  an  interpretation  of  the  importance  of  the 
zeros  throughout  the  Fresnel  region  [Fiddy  etal,  (1982)]. 

While  the  early  work  on  phase  retrieval  from  the  Fourier  transform 
modulus  data  concentrated  on  one-dimensional  functions,  the  problem 
for  the  case  of  two-dimensional  data  was  introduced  by  Napier  and  Bates 
[1974].  For  the  most  common  case  of  sampled  data,  the  problem  can  be 
stated  in  terms  of  the  ability  to  factorize  polynomials  [Bruck  and  Sodin, 
(1979)].  By  the  fundamental  theorem  of  algebra,  any  polynomial  of 
degree  two  or  more  of  a  single  variable  can  always  by  factored,  but  there 
is  no  guarantee  of  the  ability  to  factorize  polynomials  of  more  than  one 
variable.  This  basic  difference  between  the  one-dimensional  and  two- 
dimensional  problems  was  one  aspect  in  the  questions  about  ambiguity 
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and  uniqueness  of  any  reconstruction  [Bates,  (1982);  Garden  and  Bates, 
(1982);  Fright  and  Bates,  (1982)].  These  questions  will  not  be  addressed 
here,  but  rather  the  reader  is  directed  to  some  of  the  extensive  literature 
on  this  subject  [Huiser  and  vanToorn,  (1980);  Fiddy  et  a/,  (1983);  Brames, 
(1985,  1986);  Hayes  and  McClellan,  (1982);  Kierdron,  (1981);  Sanz  and 
Huang,  (1983);  Scivierand  Fiddy,  (1985)]. 

Along  with  the  many  analytic  techniques  mentioned  above,  computer 
algorithms  were  also  developed.  Most  of  the  computer  algorithms  were 
based  upon  an  iterative  procedure  first  described  by  Gerchberg  and 
Saxton  [1971,  1973].  The  original  procedure  outlined  by  Gerchberg  and 
Saxton  used  the  intensity  in  both  the  image  plane  and  the  Fourier  plane  to 
find  the  phase  function  in  each  domain.  The  algorithm  is  iterative,  and 
transforms  back  and  forth  between  the  two  domains.  At  each  iteration 
the  known  magnitude  information  is  applied  in  the  appropriate  domain, 
while  the  phase  functions  are  allowed  to  iterate  to  a  value  which  is 
consistent  with  the  known  information. 

It  should  be  noted  that  this  procedure  can  be  easily  generalized  to 
many  other  situations.  The  basic  idea  of  the  algorithm  is  simply  to 
transform  between  two  domains,  and  apply  any  known  information  in 
each  domain  as  appropriate.  One  such  generalization  was  developed  by 
Misell  [1972, 1973a,  1973b,  1973c]  for  the  case  of  information  recorded  in 
the  image  plane  and  in  one  defocused  plane.  Although  a  non-iterative 
solution  to  the  phase  problem  is  possible  when  data  from  two  different 
locations  is  known  [Dallas,  (1975);  (1976)],  Boucher  [1980a,  1980b]  has 
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shown  that  these  techniques  are  much  more  sensitive  to  noise  and  errors 
than  are  the  iterative  techniques. 

Modified  versions  of  the  Gerchberg-Saxton  algorithm  have  been 
developed  by  Fienup  [1978,  1979,  1981a,  1981b,  1982,  1983].  The  basic 
modification  to  the  algorithm  is  in  the  use  of  only  one  intensity 
measurement  (in  the  Fourier  plane),  and  the  use  of  non-negativity  and 
spatial  extent  constraints  in  the  object  domain.  Other  modifications 
involve  variations  on  how  the  known  information  and  constraints  are 
applied  in  each  domain.  Variations  of  these  reconstruction  algorithms 
require  from  60  to  200  iterations  for  a  satisfactory  solution.  A  nice 
summary  and  comparison  of  these  different  techniques  has  been  given  by 
Fienup  [1982].  Current  forms  of  Fienup's  algorithm  are  concerned  with 
overcoming  the  common  problems  of  stagnation,  improving  the 
convergence  time  (i.e.  number  of  iterations)  of  the  reconstruction  [Fienup 
and  Wackerman,  (1986);  Fienup,  (1986)],  and  working  with  complex 
valued  objects  [Fienup,  (1986a)]. 

All  of  the  reconstruction  techniques  discussed  above  have  dealt  with 
the  problem  of  phase  retrieval.  The  corresponding  problem  of  magnitude 
retrieval  was  first  investigated  by  Hayes,  Urn  and  Oppenheim,  [1980]  and 
others  [Oppenheim  et  al,  (1980);  Hayes,  (1981,  1982);  Oppenheim  et  al, 
(1982);  Oppenheim  etal,  (1983);  Oppenheim  and  Lim,  (1983)].  Although 
in  many  cases  of  interest  a  closed  form  solution  to  the  problem  of 
reconstructing  an  image  from  Fourier  transform  phase  information  is 
possible,  the  large  set  of  equations  which  must  be  solved  simultaneously 
makes  these  methods  of  solution  very  difficult.  Using  only  the  phase 
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information  of  the  Fourier  transform  as  the  known  data  in  a  Gerchberg- 
Saxton  type  algorithm  leads  to  a  recognizable  image  in  just  one  iteration, 
and  a  fairly  good  reconstruction  in  just  several  iterations  [Hayes,  (1982)]. 
Good  reconstructions  were  found  even  when  only  one  bit  of  the  phase 
information  was  combined  with  the  magnitude  information  [Van  Hove, 
(1982);  Van  Hove  etal,  (1983)]. 

The  theory  and  method  of  projections  onto  convex  sets  [Youla  and 
Webb,  (1982);  Sezan  and  Stark,  (1982)]  present  a  rigorous  mathematical 
basis  for  many  of  the  iterative  reconstruction  techniques  discussed  above. 
The  method  of  alternating  projections  onto  convex  sets  has  been  applied 
to  both  the  phase  retrieval  [Levi  and  Stark,  (1984)]  and  the  magnitude 
retrieval  [Levi  and  Stark,  (1983)]  problems. 

Other  methods  of  image  reconstruction  have  been  investigated 
[Frieden,  (1972);  Kiedron,  (1980);  Arsenault  and  Chalasinka-Macukow, 
(1983);  Bortz,  (1983);  Deighton  etal,  (1985);  Harque  and  Meyer,  (1986); 
Fienup,  (1986b)]  but  will  not  be  described  here  since  our  main  interest  is  in 
a  reconstruction  method  based  upon  the  Gerchberg-Saxton  algorithm. 

1 .3  Overview  of  this  Dissertation 

The  problem  of  reconstructing  an  unknown  image  from  partial  data  in 
the  Fresnel  zone  region  is  discussed.  The  Fresnel  region  of  a  coherent 
optical  processing  system  is  used  as  a  basis  for  the  definition  of  a  Fresnel 
zone  transform  pair.  An  Iterative  method  of  reconstructing  an  unknown 
image  located  in  the  object  domain,  from  either  the  phase-only  or 
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magnitude-only  data  in  the  transform  domain,  is  investigated  through  the 
use  of  digital  calculations. 

In  Chapt.  2,  the  Fresnel  zone  transform  is  implemented  digitally  using  a 
fast  Fourier  transform  routine.  The  number  of  Fresnel  zones,  which  is  a 
measure  of  how  far  out  of  the  Fourier  transform  plane  the  original  (or 
available)  data  is  obtained,  is  introduced  as  a  dimensionless  parameter 
used  throughout  this  dissertation.  The  number  of  Fresnel  zones  is  limited 
by  certain  sampling  considerations  which  are  discussed.  From  the  Fourier 
transform  plane  through  the  Fresnel  .  •.  'f  demonstrated  that  the 
phase-only  information  'n  the  transform  allows  for  a  very  good 
reconstruction  in  just  a  few  iterations  of  the  reconstruction  algorithm. 
While  reconstructing  an  image  from  the  magnitude  of  the  Fourier 
transform  is  very  difficult,  as  one  moves  further  into  the  Fresnel  zone 
region  good  reconstructions  are  obtained  in  only  several  iterations 
[Rolleston  and  George,  (1986)]. 

Chapter  3  presents  an  explanation  for  the  success  of  the  magnitude- 
only  reconstructions  from  data  far  into  the  Fresnel  region.  An 
approximation  of  the  Fresnel  zone  transform  integral  is  derived  which  is 
valid  for  a  large  number  of  Fresnel  zones.  The  stationary  phase 
approximations  which  are  used  are  developed  for  the  case  of  one¬ 
dimensional  functions,  with  a  summary  given  for  the  case  of  two- 
dimensional  images.  Critical  points  of  the  second  kind  are  also  calculated 
for  the  case  of  two  dimensional  images.  These  approximations  are  used 
to  show  that  the  first  iterative  guess  is  approximately  equal  to  the 
unknown  image  [Rolleston  and  George,  (1987)]. 
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As  evidenced  by  the  reconstruction  experiments  in  Chapt.  2,  phase- 
only  information  in  the  Fresnel  zone  transform  leads  to  very  good 
reconstructions  throughout  the  Fresnel  zone  region.  It  has  been  shown 
that,  for  a  given  number  of  iterations,  there  is  a  slight  degradation  in  the 
quality  of  the  reconstructed  image  as  the  number  of  Fresnel  zones  is 
increased,  but  that  the  phase  of  the  Fresnel  zone  transform  preserves  the 
structure  of  the  sharp  edges  in  the  image.  Through  several  illustrative 
examples  presented  in  Chapt.  4,  the  ability  of  the  phase  information  to 
preserve  the  edge  information  is  demonstrated  analytically. 

Chapter  5  concerns  itself  with  a  discussion  of  the  effects  of  noise  on  the 
reconstruction  process.  A  degradation  in  the  reconstructed  image  results 
if  noise  is  present  in  the  given  data.  A  model  for  coherently  additive 
noise,  and  the  resulting  probability  density  functions  of  the  magnitude 
and  phase  of  a  signal,  are  derived.  The  effects  of  both  signal  independent 
additive  and  signal  dependent  multiplicative  noise  are  investigated. 

Another  possible  source  of  error  is  the  incorrect  number  of  Fresnel 
zones  used  in  the  reconstruction  algorithm.  Chapter  6  discusses  the 
effects  of  using  the  incorrect  number  of  Fresnel  zones  and  describes  a 
method  of  searching  through  the  Fresnel  region  to  find  the  correct 
location  of  the  Fresnel  zone  transform. 

An  alternate  point  of  view  of  the  Fresnel  zone  transform,  as  defined  in 
the  context  of  the  coherent  optical  processor  discussed  in  Chapt.  2,  is 
described  in  App.  A.  This  alternate  point  of  view  allows  for  a 
generalization  of  the  Fresnel  zone  transform  to  include  the  case  of 
illumination  by  a  general  aberrated  wavefront.  An  example  of  this  type 
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Chapter  Two 


Image  Reconstruction  from  Partial  Fresnel  Zone  Information 

2.1  Introduction 

Many  of  the  methods  used  for  reconstructing  an  object  from  partial 
Fourier  transform  information  discussed  in  Chapt.  1  use  a  modified 
Gerchberg-Saxton  algorithm.  There  are  instances  when  only  the 
undistorted  phase  or  magnitude  of  the  Fourier  transform  is  known,  and  a 
reconstruction  of  the  original  object  is  desired.  It  is  generally  assumed  in 
these  reconstructions  that  the  original  object  is  real  valued. 

In  this  chapter  a  method  for  the  reconstruction  of  an  image  from 
partial  data  in  the  Fresnel  zone  region  is  presented.  The  ability  to  obtain 
excellent  image  reconstructions  using  either  phase^only  or  magnitude- 
only  data  from  the  Fresnel  zone  transform  of  a  real-valued,  non-negative 
object  will  be  demonstrated  with  the  use  of  computer «  mulations. 

The  standard  4-F  canonical  optical  processor  discussed  in  Sect.  2.2  is 
used  as  the  basis  for  the  derivation  of  the  optical  Fresnel  zone  transform 
[see  Goodman,  (1968)  §5.2  and  §7.4  for  a  further  discussion  of  coherent 
optical  processors].  From  this  groundwork  in  coherent  optical  processors, 
a  generalized  Fourier  transform  pair  is  defined.  The  numerical  methods 
used  to  calculate  the  Fresnel  zone  transforms,  as  well  as  the  sampling 
considerations,  are  discussed.  A  dimensionless  parameter,  Zp  (the  number 
of  Fresnel  zones),  is  introduced  as  a  measure  of  how  far  out  of  the  Fourier 
transform  plane  the  Fresnel  zone  transform  is  located. 
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An  alternate  point  of  view  of  the  Fresnel  zone  transform  is  described  in 
App.  A.  The  presentation  in  App.  A  formulates  the  reconstruction 
problem  in  the  context  of  Fourier  transforms.  The  problem  is  to 
reconstruct  the  magnitude  of  an  image;  given  either  the  magnitude-only 
or  phase-only  data  of  the  Fourier  transform,  and  the  phase  part  of  the 
original  complex  image.  In  this  context,  the  Fresnel  zone  transform 
discussed  in  this  dissertation  is  a  special  case  of  the  reconstruction  problem 
presented  in  App.  A. 

Section  2.3  presents  an  algorithm  used  for  reconstructing  an  image 
from  phase-only  data.  The  results  of  reconstructions  from  several 
different  Fresnel  zone  locations,  as  well  as  measures  of  the  error  in  these 
reconstructions,  are  presented.  Magnitude-only  reconstructions  are 
discussed  in  Sect.  2.4,  which  includes  error  measurements,  and  examples 
of  several  reconstructions. 

It  is  shown  that  for  this  type  of  reconstruction,  as  one  moves  away  from 
the  Fourier  transform  plane,  very  good  reconstructions  are  obtained  in 
just  a  few  iterations  using  phase-only  data.  For  the  same  number  of 
iterations,  reconstructions  from  magnitude-only  information  improve  as 
one  moves  farther  from  the  Fourier  plane. 
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2.2.  Optical  Fresnel  Zone  Transform  and  Digital  Calculation 

Consider  the  standard  4-F  canonical  optical  processor  shown  in  Fig.  2- 
1.  It  consists  of  an  input  plane  (x,y),  a  generalized  transform  plane  (u,v),  a 
Fourier  transform  plane  located  midway  between  the  two  lenses,  and  an 
output  plane  (r,s).  Throughout  this  section  we  adopt  the  notation  that 
g(x,y)  is  the  scalar  component  of  the  field  in  the  input  plane  and  G(u,v}  is 
the  scalar  component  in  the  transform  plane.  Note  that  the  transform 
plane  is  not  constrained  to  a  location  midway  between  the  two  lenses,  so 
it  is  not  necessarily  a  Fourier  transform.  If  the  input,  with  amplitude 
transmittance  g(x,y),  is  illuminated  by  a  monochromatic  plane  wave,  then 
the  scalar  component  of  the  field  in  the  transform  plane  a  distance  di 
behind  the  first  lens  will  be  given  by 


f  e 

G(u,i>)  =  Cj  g(.x^)K^(xj^-,u,v)dxdy  , 

J  J  —00 

where  the  transfer  kernel  Ki(x,y;  u,v)  is  as  follows: 

-ina  <*^+y^)  -i2B{f  x  *  f  y) 

K^(xjr;  u.v)  =  e  e  ^  . 

For  simplicity  of  notation,  we  have  defined  the  complex  constant  ci 


(2.2-1) 


(2.2-2) 


by 


(« 


-  i*(dj  *-Fi 


"l  = 


\F 


and  the  spatial  frequency  variables  fx  and  fy  by 

f  =  —ul\F  and  f  =  —vl\F  . 

'  X  y 


(2.2-3) 


(2.2-4) 


The  convenient  offset  parameter  ai,  used  in  Eq.  (2.2-2),  is  defined  as 
follows; 


- 17- 


(2.2-5) 


ifdi  sF,  thenai  sO.  and  the  well  known  Fourier  transform  relationship 
between  the  front  and  back  focal  planes  of  a  convex  lens  is  obtained.  If 
di=^  F,  then  the  generalized  transform  plane  is  in  the  Fresnel  zone  of  the 
optical  processor;  this  is  the  region  of  central  interest  in  this  discussion. 
The  terminology  optical  Fresnel  zone  transform  will  be  used  to  describe 
Eq.  (2.2-1). 

If  the  scalar  component  of  the  field,  G(u,v),  is  propagated  a  distance  d2 
to  the  second  lens  and  then  to  the  output  plane  (r,s),  the  scalar 
component  of  the  output  is  given  by 


II  /'■>  G(u,v) K^(u,v;  r,s)  du  dv  , 

where  the  transfer  kernel  K2(u,v;  r,s)  Is  as  follows: 

-ina,(r*+»^)  -i2n(/  r*f  «) 
K2(u,u;r,8)  =  e  ^  e  *  ^  . 

The  complex  constant  C2  is  given  by 


(2.2-6) 


(2.2-7) 


-ikid  (2.2-8) 

le 

XF  ’ 

and  the  spatial  frequencies  fx  and  fy  are  defined  in  Eq.  (2.2-4).  A 
convenient  offset  parameter  02.  used  in  Eq.  (2.2-6),  Is  defined  as  follows: 


For  the  case  when  di  +  d2  =  2F,  the  output  gout(r,s)  is  given  by 


(2.2-9) 


^^(r,s)=-c 


(2.2-10) 
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Aside  from  a  phase  retardation  term,  the  output  is  an  inverted  image  of 
the  input. 


2.2.2  The  Fresnel  Zone  Transform  Pair 

it  is  convenient  to  summarize  the  above  discussion  with  the  following 
definition  of  a  Fresnel-zone  transform  pair.  Consider  a  general  function, 
g(x,y),  and  its  associated  Fresnel  zone  transform  G(fx,  fy).  The  Fresnel-zone 
transform  of  g(x,y)  is  given  by 


f  ,f  )dxdy  , 

OD  ^ 


and  the  inverse  Fresnel  zone  transform  is  given  by 


(2.2-11) 


f  +®  f  +• 

g(x.y)=  I 

The  Fresnel-zone  transform  kernel  K  (x,y;  fx,fy)  is  given  as 

2  2 

>iaa,(s  %'¥f  y) 

Kix.y;f^.f^)  =  e  ‘  e  *  ^  . 


(2.2-12) 


(2.2-13) 


and  *  denotes  a  conjugate  complex.  In  this  mathematical  definition  of  a 
Fresnel-zone  transform  pair  the  phase  retardation  and  image  inversion 
that  occur  in  the  optical  transform  have  been  suppressed. 

This  Fresnel-zone  transform  pair  is  very  closely  related  to  the  well 
known  Fourier  transform  pair.  A  Fresnel  zone  transform  is  easily 
calculated  by  using  the  analytic  properties  of  Fourier  transform  theory  as 
well  as  the  numerical  methods  associated  with  the  fast  Fourier  transform. 
To  calculate  the  Fresnel  zone  transform  of  an  input  object  g(x,y),  first 
multiply  the  object  by  the  complex  Gaussian  phase  term,  expl-inai(x2  + 
y2)],  and  then  take  the  Fourier  transform  of  this  product.  This 


multiplication  followed  by  a  Fourier  transform  is  much  easier  and  faster 
than  convolving  the  Fourier  transform  of  the  object  with  the  Fourier 
transform  of  the  complex  Gaussian  phase  term.  The  Fresnel  zone 
transform  of  a  digitized  image  can  be  calculated  easily  using  a  fast  Fourier 
transform  (FFT)  routine. 

If  the  Fresnel  zone  transform  G(fx,fy)  is  known,  then  an  inverse  Fresnel 
zone  transform  can  be  used  to  find  the  unknown  image.  To  calculate  an 
inverse  Fresnel  zone  transform,  first  perform  an  inverse  Fourier 
transform,  and  then  multiply  by  the  complex  Gaussian  phase  term 
exp[inai(x2  -i-  y2)].  If  it  is  known  that  the  object  is  real  valued  and  non¬ 
negative,  then  the  modulus  of  the  inverse  Fourier  transform  can  be  used 
as  the  inverse  Fresnel-zone  transform  rather  than  multiplying  by  the 
complex  Gaussian.  The  inverse  Fresnel  zone  transform  can  be  calculated 
using  an  inverse  fast  Fourier  transform  (IFFT)  routine,  followed  by  a 
multiplication  or  absolute  value  operation. 

Throughout  this  dissertation  the  Fresnel  zone  transform  pair  defined 
in  Eqs.  (2.2-11)  and  (2.2-12)  will  be  used,  rather  than  the  optical  transform 
described  in  Eqs.  (2.2-1)  and  (2.2-6).  This  prevents  the  resultant  image 
from  being  inverted  and  suppresses  the  constant  phase  retardation  that 
occurs  in  the  optical  transform.  Since  the  input  functions  are  limited  to 
real  valued  non-negative  objects,  the  complex  Gaussian,  exp[inai(x2  +  y2)], 
used  to  calculate  an  inverse  Fresnel  zone  transform  is  not  needed. 
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2.2.3  Sampling  Considerations 

In  the  numerical  computation  of  the  Fresnel-zone  transform  defined  in 
Eq.  (2.2-1 1)  a  consideration  of  the  sampling  interval  must  be  made.  The 
sampling  interval  is  controlled  both  by  the  spectrum  of  g(x,y)  and  by  the 
oscillatory  nature  of  the  complex  Gaussian  multiplier.  This  discussion  will 
investigate  the  sampling  limits  based  upon  the  oscillations  of  the  complex 
Gaussian  term,  and  not  on  the  spectrum  of  g(x,y). 

Consider  the  s^m  ilmg  of  the  input  in  more  detail.  From  Eqs.  (2.2-11) 
and  (2.2-13),  it  ;;  seen  that  the  process  of  calculating  a  Fresnel  zone 
transform  is  the  same  as  calculating  the  Fourier  transform  of  the  product 

-in^AxUy^)  (2.2-14) 

g(x,y)  e 

If  this  function  is  sampled  at  xsmT  and  y  =  nT,  where  T  is  the  sampling 
interval,  ms:-N/2.  -N/2  +  1, ... ,  +  N/2-1,  and  n  =-N/2, -N/2  +  1, ....  +  N/2-1, 
Eq.  (2.2-14)  becomes 


-I  n  a,  T^(  m  ^  +  n  ^ )  (2.2-15) 

g(mT,  nT)  e 

The  total  number  of  samples  is  NxN.  An  FFT  can  be  used  on  this  sampled 
data  only  if  the  function  has  been  sampled  at  a  sufficient  frequency.  The 
term  g(mT,  nT)  will  be  a  digitized  picture  of  finite  extent  and  it  is  assumed 
that  T  can  be  chosen  so  that  the  finest  detail  of  interest  in  the  image  will 
be  properly  sampled.  However,  the  term  exp[-inaiT2(m2  +  n2)]  can 
oscillate  very  rapidly  near  the  edge  of  the  image.  To  obtain  an  estimate  of 
how  large  ai  can  be,  let  n=0  and  assume  that  the  argument  of  the 
complex  Gaussian  is  to  change  by  no  more  than  n/4  from  one  sample  point 
to  the  next  at  the  edge  of  the  image.  Under  these  constraints,  the 
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following  inequality  nriust  be  satisfied: 


J  -N 

21 

-N  \2 

n 

na,TT  - 

)  -  n  q.T* 

-  +  1  ) 

£  - 

1  \  2 

/  > 

1  2  / 

4 

(2.2-16) 


With  N  >  >  1  this  reduces  to 


I  aj  s 


1 

4NT^ 


(2.2-17) 


Combining  Eqs.  (2.2-5)  and  (2.2-17)  we  find  that  the  distance  di  must 
satisfy  the  constraint; 


\F 

s  - . 

4NT* 


(2.2-18) 


These  constraints  are  based  upon  the  assumption  that  the  argument  of 
the  complex  Gaussian  should  change  by  no  more  than  n/4  from  one 
sample  point  to  the  next.  To  compare  this  sampling  criterion  with  the 
Nyquist  sampling  rate,  consider  the  real  part  of  the  complex  Gaussian 
given  in  Eq.  (2.2-15): 

cos(naj  T^m^),  (2.2-19) 

where  along  the  x-axis,  n  is  equal  to  zero.  This  expression  can  be  written 
as 


cos(2n  ( -^  T  ^  ml  m) , 

ft 


(2.2-20) 


where  the  term  in  square  brackets  can  be  viewed  as  an  instantaneous 
frequency,  fi. 

At  the  edge  of  the  image  (i.e.  n  =  ±  N/2),  the  instantaneous  frequency 
is  at  a  maximum  and  is  given  by 


(2.2-21) 


The  Nyquist  sampling  interval  of  a  sinusoid  of  this  frequency  is  at  least 
1/(2fi),  but  since  the  sampling  interval  is  already  given  by  Am  =  1,  there 
will  be  a  limit  on  the  magnitude  of  fi.  The  Nyquist  limit  on  this 
instantaneous  frequency  is  given  by 


^  2 

Combining  Eqs.  (2.2-21)  and  (2.2-22)  gives  an  upper  bound  onai, 

I  .  2 

lajs  - 

*  A 


(2.2-22) 


(2.2-23) 


which  is  eight  times  the  limit  given  by  Eq,  (2.2-17).  This  factor  of  eight  is 
the  result  of  using  the  change  in  the  total  argument  of  the  complex 
Gaussian,  and  limiting  this  change  to  n/4  rather  then  n  in  the  derivation  of 
Eq.  (2.2-17). 


2.2.4  Definition  of  Zf 

In  going  from  the  center  to  the  edge  of  the  image,  the  argument  of 
the  complex  Gaussian  in  Eq.  (2.2-15)  either  increases  or  decreases 
quadratically,  depending  upon  the  sign  of  ai .  Each  multiple  of  2n  that  the 
argument  goes  through  will  be  counted  as  one  Fresnel  zone.  Along 
either  the  x  or  y  axis  the  maximum  value  of  this  argument  will  be 


.1., !.■(?)' 


(2.2-24) 


Dividing  this  quantity  by  2n  gives  the  total  number  of  Fresnel  zones.  For 
any  given  value  of  at ,  as  defined  in  Eq.  (2.2-5),  the  correspo  <ding  number 
of  Fresnel  zones,  Zf,  is  given  by; 
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“,/NT\2  (2.2-25) 

The  dimensionless  parameter  Zf  will  be  used  throughout  this  dissertation 
to  denote  how  far  into  the  Fresnel  zone  the  transform  plane  is  located.  It 
is  important  to  remember  that  as  the  magnitude  of  Zf  increases,  the 
transform  plane  is  moving  farther  away  from  the  Fourier  transform  plane, 
and  further  into  the  Fresnel  zone  region  of  the  optical  processor. 

Using  the  maximum  value  for  ai  given  in  Eq.  (2.2-17),  the  argument  of 
the  complex  Gaussian  changes  from  zero  to  2n(N/32)  in  moving  from  the 
center  to  the  edge  of  the  image.  The  maximum  number  of  Fresnel  zones 
that  can  be  sampled  is  thus  defined  as  N/32.  The  constraint  on  di  given  in 
Eq.  (2.2-18)  is  equivalent  to 
Fresnel  zones,  Zf,  satisfy  the 


a  constraint  requiring  that  the  number  of 


equation: 
N 


32 


(2.2-20) 
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25 


2.3.  Reconstruction  from  Fresnel  Zone  Phase  Information 

An  iterative  method  of  reconstructing  an  object  from  Fresnel  zone 
transform  phase  information  will  be  discussed  in  this  section.  In  this 
reconstruction  technique  it  is  assumed  that  the  phase  of  the  Fresnel  zone 
transform  is  known  and  that  the  original  object  is  real  and  non-negative. 
It  is  also  assumed  that  the  phase  information  is  free  from  noise  and  that 
the  value  of  Zf  is  known,  i.e.  the  location  in  the  Fresnel  zone  region  at 
which  the  phase  information  was  obtained  is  known.  A  block  diagram 
illustrating  this  iterative  algorithm  is  shown  in  Fig.  2-2. 

The  algorithm  starts  with  an  NxN  object  of  random  values,  taken  from 
a  uniform  distribution  with  a  range  of  0  to  255,  as  the  initial  guess  for  the 
object.  The  guess  for  the  object  is  multiplied  by  the  known  complex 
Gaussian  phase  term  and  then  padded  with  zeros  to  a  size  of  2Nx2N.  An 
FFT  is  performed  to  calculate  a  guess  for  the  Fresnel-zone  transform.  The 
modulus  of  this  guess  is  multiplied  by  the  known  phase  information  and 
an  IFFT  is  performed.  Taking  the  modulus  of  the  central  NxN  pixels  yields 
the  new  guess  for  the  object.  The  algorithm  continues  by  transforming 
back  and  forth  between  the  object  domain  and  the  Fresnel  zone  while 
applying  the  known  information  or  constraints  in  each  domain. 

A  direct  measure  of  the  error  between  the  known  phase  information 
and  the  phase  of  the  Fresnel  zone  transform  of  the  current  guess  is  given 
by  the  absolute  phase  error,  Ep(*).  The  absolute  phase  error  is  defined  by 


N-l  N-l 
1  1 


m=0  n=0 


exp{i  d(m,n)}  —  exp{i  0^  (m,n)} 


(2.3-1) 
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where  exp{i6(m,n)}  is  the  known  phase  information,  and  exp{i6i(m,n)}  is 
the  phase  of  the  Fresnel  zone  transform  of  the  guess  for  the  image  after  i 
iterations. 

At  each  iteration  an  analytic  measure  of  the  relative  squared  error  in 
the  image  plane,  Ei(i),  is  given  by 


S-l  N-l 

1  1 

m^O  n=0 


g  Jim, n)  -  g^im.n) 


N-l  N-l 

1  1 


gim.n) 


(2.3-2) 


in*0  »-0  ’ 

where  i  is  the  iteration  number,  gi(m,n)  is  the  calculated  guess  for  the 
image  after  i  iterations  and  before  the  image  plane  constraints  have  been 
imposed,  and  gj'(m,n)  is  the  result  of  applying  the  image  plane  constraints 
on  gi(m,n).  Although  this  error  is  very  difficult  to  relate  to  the  known 
transform  magnitude  information,  this  error  does  give  an  indication  as  to 
how  much  the  guess  for  the  image  must  be  modified  at  each  iteration  to 
fulfill  the  image  plane  constraints. 

Figure  2-3  shows  the  original  256x256x8-bit  digitized  image  used  in  the 
reconstruction  experiments.  The  structure  of  the  two  railroad  cars  on  the 
left-hand  side  of  the  picture  provide  a  good  reference  for  two  different 
spatial  frequencies,  ls  does  the  low  contrast  fine  detail  at  the  back  of  the 
railroad  car  in  the  center  of  the  picture.  The  uniform  illumination  of  the 
sky  and  the  hard  edge  of  the  coal  pile  against  the  background  of  the  sky 
are  also  good  reference  structures. 

A  256x256  complex  Gaussian  was  calculated  for  several  different 
Fresnel  zone  transforms.  The  known  phase  information  that  was  used  in 


-27- 


Fig.  2-3  Original  digitized  image  used  in  the  reconstruction 
experiments.  The  image  has  256  X  256  pixels  with  8 
bits/pixel. 

the  iterative  reconstruction  vy^as  calculated  by  multiplying  the  known 
image  by  the  corresponding  complex  Gaussian,  padding  this  product  with 
zeros  to  obtain  a  512x512  image,  performing  an  FFT  with  this  data,  and 
then  setting  the  magnitude  equal  to  unity  at  every  pixel.  This  known 
Fresnel  zone  phase  information  from  four  different  Fresnel  zone 
transforms  is  shown  in  Fig.  2-4.  This  phase  information  was  calculated  for 
Zf  =  0.0.  0.5.  3.0,  and  8.0  and  is  shown  in  Figs.  2-4(a).  2-4(b).  2-4(c).  and  2- 
4(d)  respectively.  The  phase  is  represented  by  256  grey  levels,  with  a  phase 
which  is  just  greater  than  -  n  corresponding  to  black,  and  a  phase  which  is 
equal  to  n  being  shown  as  white. 

The  reconstructions  for  the  object  after  the  fifth  iteration  using 
different  values  of  di  are  shown  in  Fig.  2-5.  The  reconstruction  in  Fig.  2- 
5{a)  corresponds  to  Zf  =  0.0.  which  is  the  Fourier  transform  plane. 


-29- 


(a)  Zf  = 


(b)ZF  =  0.5 


0.0 


(c)  Zf  =  3.0  (d)  Zf  =  8.0 

Fig.  2-4  Known  phase  information  from  four  different  Fresnel  zone  transforms.  As 
defined  in  Eqs.  (2  2-5)  and  (2.2-25),  Zf  is  the  number  of  Fresnel  zones,  or  a 
measure  of  the  distance  away  from  the  Fourier  plane.  Phase  values  from  -  n  to 
I  +  n  are  mapped  to  a  grey  scale  from  black  to  white. 


Reconstructions  limited  to  the  Fourier  transform  plane  have  been 
discussed  in  prior  work  by  Hayes,  Lim,  and  Oppenheim  [1980,  1982]  and 
Hayes  [1982],  where  they  demonstrated  that  a  good  reconstruction  is 
obtained  after  only  one  iteration.  It  has  been  found  that  good 
reconstructions  are  obtained  in  only  one  iteration  from  phase-only  Fresnel 
zone  data  as  well.  For  the  purpose  of  comparing  phase-only 
reconstructions  with  magnitude-only  reconstructions,  the  results  of 
reconstructions  after  five  iterations  are  shown. 

Figures  2-5(b),  2-5(c),  and  2-5(d)  show  the  reconstructions  in  the  Fresnel 
zone  region.  These  reconstructions  are  shown  with  Zf  =  0.5,  Zf  =  3.0,  and 
Zf  =  8.0  respectively.  Note  that  an  increase  in  the  number  of  Fresnel  zones 
corresponds  to  an  increase  in  the  distance  away  from  the  Fourier 
transform  plane.  The  high  frequency  detail  and  edges  are  very  well 
preserved  in  all  of  these  reconstructions.  As  the  number  of  Fresnel  zones 
increases,  the  uniformity  of  the  sky  is  degraded,  but  overall  there  is  still  a 
very  good  reconstruction.  From  these  results  we  see  that  this  type  of 
reconstruction  is  fairly  insensitive  to  a  shift  out  of  the  Fourier  transform 
plane. 

Figure  2-6  is  a  plot  of  the  absolute  phase  error,  Ep(>),  for  iterations  2 
through  20  at  several  different  Fresnel  zone  locations.  The  error  at  the 
first  iteration  is  approximately  .5,  and  is  omitted  here  to  provide  for  a 
more  detailed  vertical  scale.  There  is  a  slight  increase  in  the  absolute 
phase  error  as  the  number  of  Fresnel  zones  increases;  this  agrees  with  the 
slight  degradation  in  the  visual  quality  of  the  images  shown  in  Fig.  2-5. 
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A  plot  of  the  relative  squared  error  in  the  image  plane  at  several 
different  Fresnel  zone  locations  is  shown  for  iterations  2  through  20  in  Fig. 
2-7.  The  error  after  the  first  iteration  is  approximately  .10  to  .14  and  is 
omitted  in  this  plot  to  allow  for  greater  resolution  in  the  vertical  scale. 
There  is  very  little  difference  in  the  error  as  the  transform  plane  moves 
away  from  the  Fourier  plane  and  further  into  the  Fresnel  zone  region. 
The  relative  squared  error  in  the  image  plane  is  independent  of  the 
number  of  Fresnel  zones. 


Fig.  2-6  Error  in  transform  plane,  Ep('),  vs.  iteration  number  for  phase-only 
reconstruction  at  Fresnel  zones  equal  to  0.0,  0.5,  3.0,  and  8.0. 
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Fig.  2-7  Error  in  image  plane,  vs.  iteration  number  for  phase-only  reconstructions 
at  Fresnel  zones  equal  to  0.0, 0.5, 3.0,  and  8.0. 

As  the  reconstruction  algorithm  progresses,  there  is  a  steady  decrease 
in  both  the  image  plane  error  and  the  transform  plane  etror.  The 
reconstruction  after  the  first  iteration  is  an  edge  enhanced  (almost  to  the 
point  of  being  binary)  outline  of  the  unknown  image.  As  the 
reconstruction  progresses,  the  low  frequency  areas  of  the  image  are 
slowly  filled  in.  In  reconstructing  an  image  from  phase-only  data,  the 
quality  of  the  reconstruction  is  essentially  independent  of  the  number  of 
Fresnel  zones.  In  locations  throughout  the  Fresnel  zone,  the  high  contrast. 
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high  frequency  edges  which  dominate  the  first  few  iterative  guesses  are 
slowly  worn  down  so  that  the  low  frequency,  low  contrast  areas  of  the 
image  may  be  recovered.  The  slight  degradation  in  the  quality  of  the 
reconstructed  image  as  the  number  of  Fresnel  zones  increases  is  reflected 
in  a  small  increase  in  the  value  of  the  absolute  phase  error. 
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2.4  Reconstruction  from  Fresnel  Zone  Magnitude  information 

The  iterative  method  of  reconstructing  an  object  from  Fresnel  zone 
magnitude-only  information  is  very  similar  to  the  algorithm  discussed 
earlier  for  reconstructing  an  object  from  Fresnel  zone  phase-only 
information.  The  only  difference  is  in  the  information  that  's  known 
about  the  Fresnel  zone  transform.  In  this  reconstruction  it  is  assumed  that 
the  magnitude  of  the  Fresnel  zone  transform  is  known,  and  that  the 
original  object  is  real  and  non-negative.  It  is  also  assumed  that  we  know 
where  in  the  Fresnel  zone  region  the  magnitude  information  was 
calculated.  A  block  diagram  showing  this  iterative  algorithm  is  presented 
in  Fig.  2-8. 

This  iteration  also  starts  with  an  initial  random  guess  for  the  object. 
The  earlier  algorithm  is  modified  so  that  for  this  reconstruction,  the 
known  magnitude  of  the  Fresnel  zone  transform  is  multiplied  by  the 
current  calculated  guess  for  the  Fresnel  zone  transform  phase.  An  IFFT  is 
performed,  and  the  modulus  of  the  central  NxN  pixels  yields  a  new  guess 
for  the  object  in  the  same  manner  as  the  previous  algorithm. 

At  each  iteration,  an  analytic  measure  of  the  relative  squared  error  in 
the  image  plane  (as  defined  by  Eq.  2.3-2)  is  calculated.  The  absolute  error 
in  the  transform  plane,  is  given  by 

^  m  =  0  n  =  0 

where  i  is  the  iteration  number,  |G(m,n)|  is  the  measured  magnitude  data, 
and  |Gi(m,n)|  is  the  magnitude  data  calculated  from  the  Fresnel  zone 
transform  of  the  guess  for  the  image  after  i  iterations.  This  error  is  a  direct 
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guess:  g(0) 


Multiply  by  complex  Gaussian 
exp  {  -  i  X  ai  (x2  ♦  y2 ) } 
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new  guess:  gp^^i 


Fig.  2-8  Block  Diagram  showing  the  iterative  technique  of  image  reconstruction  from 
Fresnel  zone  magnitude-only  information. 
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measure  of  how  well  the  iteration  is  progressing  in  terms  of  being 
consistent  with  the  known  magnitude  information. 

The  original  image  shown  in  Fig.  2-3  was  also  used  in  this  set  of 
experiments.  The  known  magnitude  information  was  calculated  by 
multiplying  the  image  by  the  correct  complex  Gaussian,  padding  the  data 
with  zeros,  calculating  the  FFT,  and  then  setting  the  phase  equal  to  zero  at 
each  pixel.  The  known  Fresnel  zone  magnitude  information  from  four 
different  Fresnel  zone  transforms  is  shown  in  Fig.  2-9.  This  magnitude 
information  was  calculated  for  Zf  =  0.0,  0.5,  3.0,  and  8.0  and  is  shown  in 
Figs.  2-9(a),  2-9(b),  2-9(c),  and  2-9(d)  respectively.  The  magnitude  is 
displayed  with  256  grey  levels  such  that  the  maximum  value  of  loge(1. 

|G| )  is  shown  as  white  and  |G|  =  0  is  shown  as  black.  The  dynamic  range  of 
the  computer  is  much  larger  than  zero  to  255,  so  this  logarithmic  scaling 
gives  a  much  better  representation  of  the  data  than  a  direct  linear  scaling 
of  the  magnitude.  From  this  figure  we  can  see  that  there  is  a  spreading  of 
the  transform  pattern  as  we  move  further  into  the  Fresnel  region. 

The  reconstructions  for  the  object  after  the  fifth  iteration  using 
different  values  of  di  are  shown  in  Fig.  2-10.  The  reconstruction  in  Fig.  2- 
10(a)  corresponds  to  Zf  =  0.0,  which  is  the  Fourier  transform  plane.  More 
sophisticated  methods  can  be  used  to  reduce  the  number  of  iterations 
required  for  a  given  quality  of  reconstruction,  but  will  not  be  investigated 
here.  In  the  results  discussed  here,  we  are  concerned  principally  with  the 
effect  of  moving  into  the  Fresnel  region,  and  not  with  the  use  of  different 
iterative  algorithms. 


(c)  Zf  =  3.0  (d)  Zf  =  8.0 

9.  Known  magnitude  information  from  four  different  Fresnel  zone  transforms. 
Data  is  displayed  such  that  the  maximum  value  of  loge(1.  +  IGI  )  is  shown  as 
white,  and  IGI  =  0  is  shown  as  black. 
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(c)Zf  =  3.0  {d)ZF  =  8.0  ^ 

Fig.  2-10.  Results  of  reconstructing  an  object  from  magnitude-only  data  at  four 
different  locations  in  the  Fresnel  zone  region.  The  reconstructions  are 
shown  after  five  iterations.  As  we  move  further  away  from  the  Fourier 

transform  plane  the  reconstructions  improve.  j 
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Figures  2-10(b).  2-10(c),  and  2-10(d)  show  the  reconstructions  as  we 
move  further  away  from  the  Fourier  plane.  These  reconstructions  from 
the  Fresnel  zone  region  are  shown  for  Zf  =  0.5,  Zf  =  3.0,  and  Zf  =  8.0 
respectively.  Recall  that  an  increase  in  the  number  of  Fresnel  zones 
corresponds  to  an  increase  in  the  quantity  |1-di/F|.  The  reconstruction  at 
0.5  Fresnel  zones  shows  the  location  of  the  high  contrast  edges  as  well  as 
the  uniform  illumination  of  the  sky.  At  3.0  Fresnel  zones,  some  of  the  high 
frequency  detail  is  starting  to  appear  and  the  high  contrast  edges  are  well 
defined.  There  is  a  very  good  reconstruction  at  Zf  =  8.0.  In  this 
reconstruction  much  of  the  high  frequency  detail,  as  well  as  the  low 
contrast  fine  detail  at  the  back  of  the  railroad  car  in  the  center  of  the 
picture,  have  been  resolved  .  From  these  results  we  see  that  this  type  of 
algorithm  gives  a  better  reconstruction  as  the  transform  plane  moves 
further  away  from  the  Fourier  transform  plane,  and  further  into  the 
Fresnel  zone  region  of  the  optical  processor. 

Figure  2-11  is  a  plot  of  the  absolute  error  in  the  transform  plane  at 
several  different  Fresnel  zone  locations  for  iterations  2  through  75.  The 
error  changes  the  most  in  the  first  several  iterations;  from  this  point  on 
the  rate  of  change  in  the  error  is  directly  proportional  to  the  number  of 
Fresnel  zones.  After  75  iterations,  the  error  at  8.0  Fresnel  zones  has 
decreased  rapidly  to  near  zero  value,  at  3.0  Fresnel  zones  the  error  is  still 
decreasing  steadily,  and  at  0.5  Fresnel  zones  the  error  is  decreasing  very 
slowly.  The  error  at  0.0  Fresnel  zones  (the  Fourier  transform  plane)  has 
stagnated  after  approximately  50  iterations. 


i 
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Fig.  2-11  Error  in  transform  plane,  Ep<'),  vs.  iteration  number  for  Fresnel  zones  equal  to 
0.0, 0.5,  3.0,  and  8.0. 

A  plot  of  the  relative  squared  error  for  iterations  2  through  75  at 
several  different  Fresnel  zones  is  shown  in  Fig.  2-12.  The  error  decreases 
the  most  during  the  first  several  iterations  at  all  Fresnel  zone  locations. 
The  value  of  the  error  after  the  first  iteration  is  inversely  proportional  to 
the  number  of  Fresnel  zones.  For  a  given  number  of  iterations,  the  error  is 
lower  for  a  larger  number  of  Fresnel  zones.  This  agrees  well  with  the 
reconstructions  shown  in  Fig.  2-10. 
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Fig,  2-12  Error  in  image  plane,  Ei^').  vs.  iteration  number  for  Fresnel  zones  equal  to  0.0,  V| 

0.5, 3.0,  and  8.0. 

The  quality  of  the  image  in  Fresnel  zone  magnitude-only 
reconstructions  is  not  only  a  function  of  the  error  in  the  image  and 
transform  domains,  but  also  depends  on  the  rate  of  change  in  the  * 

absolute  error  in  the  transform  domain. 


Chapter  Three 


Stationary  Phase  Approximations  in  Fresnel  Zone 
Magnitude-only  Reconstructions 


3.1  Introduction 

tt  was  shown  in  Chapt.  2  that  there  is  a  very  noticeable  improvement  in 
the  quality  of  the  magnitude-only  reconstructions  as  the  number  of 
Fresnel  zones  increases.  Moving  a  small  distance  into  the  Fresnel  region 
aids  in  the  reconstruction  process  in  that  the  number  of  iterations  needed 
to  produce  a  good  image  is  fewer  than  in  the  case  of  reconstructing  an 
image  from  Fourier  transform  magnitude-only  data.  Data  from  deep  in 
the  Fresnel  zone  region  leads  to  good  reconstructions  in  just  a  few 
iterations. 

An  analytic  explanation  of  why  data  obtained  with  a  large  number  of 
Fresnel  zones  yields  a  good  reconstruction  in  just  a  few  iterations  is  the 
subject  of  this  chapter.  The  analytic  techniques  used  in  the 
approximations  discussed  below  are  part  of  the  general  area  of 
asymptotic  expansions.  A  good  general  discussion  of  asymptotic 
expansions  is  given  by  Erdelyi  [1956]  and  by  Copson  [1965],  while  a  more 
applicable  discussion  concerning  the  case  of  two-dimensional  integrals  is 
given  by  Born  and  Wolf  [1980].  The  case  of  two-dimensional  diffraction 
integrals  has  also  been  discussed  by  VanKampen  [1949, 1950]  and  Seligson 
[1981], 
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In  Sect.  3.2  the  definition  of  the  Fresnel  zone  transform  is  formulated 
so  that  the  integral  can  be  approximated  by  the  method  of  stationary 
phase.  This  method  of  approximating  the  Fresnel  zone  transform  is  then 
used  to  form  an  estimate  of  the  first  guess  for  the  Fresnel  zone  transform. 
The  first  iterative  guess  for  the  unknown  image  is  shown  to  be 
approximately  equal  to  the  unknown  image. 

The  discussion  in  Sect.  3.2  uses  the  one-dimensional  Fresnel  zone 
transform.  A  summary  for  the  case  of  two-dimensional  images  is  given  in 
Sect.  3.3.  While  the  analytic  expressions  given  in  Sects.  3.2  and  3.3  lead  to 
a  better  understanding  of  the  process  of  reconstructing  an  image  from 
Fresnel  zone  magnitude-only  data.  Sect.  3.4  presents  a  geometrical 
derivation  of  the  stationary  phase  point,  and  discusses  the  existence  of  a 
stationary  phase  point. 

Critical  points  of  the  second  kind  are  calculated  in  Sect.  3.5,  and  the 
associated  term  in  the  asymptotic  expansion  of  the  Fresnel  zone  transform 
integral  is  derived.  The  first  term  in  the  asymptotic  expansion,  due  to  the 
critical  points  of  the  first  kind  discussed  in  Sects.  3.2  and  3.3,  are  of  the 
order  (Zf)-i,  while  the  term  in  the  asymptotic  expansion  due  to  critical 
points  of  the  second  kind  goes  as  (2f)-3/2. 

An  illustrative  example  of  the  Fresnel  zone  transform  for  a  large 
number  of  Fresnel  zones,  and  the  magnitude-only  reconstruction  after 
one  iteration  is  given  in  Sect.  3.6.  The  structure  of  the  Fresnel  zone 
transform  is  shown  to  be  dominated  by  the  effects  of  critical  points  of  the 
first  kind. 
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3.2.  Critical  Points  of  the  First  Kind 

In  this  section  only  the  one-dimensional  case  of  the  Fresnel  zone 
transform  will  be  considered.  The  results  and  derivations  presented  in  the 
following  discussion  are  easily  extended  to  two  dimensions,  and  a 
summary  of  the  results  for  the  two-dimensional  case  are  given  in  Sect.  3.3. 

Consider  a  general  function,  g(x),  and  its  associated  Fresnel  zone 
transform  G(fx).  The  Fresnel  zone  transform  of  g(x)  is  given  by 

(3  2-1) 

G(f^)=  J  g(x)K{x.  f^)dx  V  i 

By  direct  calculation,  one  can  readily  verify  the  following  inversion 
relationship: 


Gif  )K\x-,f  )df  <3.2-2 

_  oe 

Where  *  denotes  a  complex  conjugate,  and  the  Fresnel-zone  transform 
kernel  K  (x;  fx)  is  given  as 


-inoJC^  -i2nfx  (3.2-3) 

Kix;f^)  =  e  c  *  . 

The  constant  a  is  related  to  the  distance  away  from  the  Fourier  plane.  The 
generalized  frequency,  fx,  will  be  equal  to  the  Fourier  spatial  frequency 
only  when  a  =  0. 

Consider  a  function  g(x)  that  is  non-zero  only  over  a  finite  region,  say 

g(x)  =  0  for  |x|  >  6  ;  (3.2-4) 
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then  the  complex  Gaussian  term  in  the  transform  kernel  will  oscillate  a 
finite  number  of  times  as  we  move  from  x  =  0  to  x  =  b.  The  argument  of 
the  complex  Gaussian  must  be  continuously  increasing  or  decreasing, 
depending  upon  the  sign  of  a.  The  number  of  multiples  of  2n  that  this 
argument  goes  through  is  the  number  of  Fresnel  zones,  Zp,  associated 
with  a  particular  Fresnel  zone  transform.  An  increase  in  the  magnitude  of 
Zp  corresponds  to  a  location  in  the  Fresnel  zone  region  which  is  farther 
away  from  the  Fourier  plane.  The  number  of  Fresnel  zones  is  given  by; 


or 


(3.2-5) 


(3.2-6) 


Combining  Eqs.  (3.2-1),  (3.2-3),  and  (3.2-6),  the  expression  for  a  Fresnel 
zone  transform  with  Zp  Fresnel  zones  can  be  expressed  in  terms  of  the 
number  of  Fresnel  zones  as; 


-in*^  2Z„/fr^  -i2nf  z 

G(/')=|  g(x)  e  e  *dx. 

'  '  -6 


(3.2-7) 


The  bounds  of  integration  reflect  the  non-zero  region  of  g(x)  defined  by 
Eq.  (3.2-4). 

To  formulate  the  expression  for  evaluation  by  the  stationary  phase 
approximation,  rewrite  Eq.  (3.2-7)  as 
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(3.2-8) 


G(/)  =  e  *  ^  g{x)e 

J  -h 


-iZ^dnt h  ) 


V  2Z^} 


dx 


3.2.2.  Stationary  Phase  Approximation 

In  this  application,  the  principle  of  stationary  phase  makes  use  of  the 
fact  that  if  Zp  is  large  enough,  and  g(x)  is  slowly  varying  compared  with 
the  complex  exponential,  then  the  value  of  the  integral  in  Eq.  (3.2-8)  is 
determined  primarily  by  the  behavior  of  the  integrand  at  certain  critical 
points.  Critical  points  of  the  first  kind,  which  will  be  considered  in  this 
section,  are  very  similar  for  both  one-  and  two-dimensional  integrals.  An 
analysis  of  critical  points  of  the  second  kind  will  be  discussed  in  Sect.  3-3.2. 

Let  the  argument  of  the  complex  exponential  be  written  as 


(3.2-9) 


where 

6V  ^ 

<|)(x)=  -s^fi[Z^(2n/62)(x-p  ,  (3.2-10) 

and  sgn[ZF]  is  equal  to  the  sign  of  Zp.  Critical  points  of  the  first  kind  are 
located  in  regions  where 


and 


U1  <  ft  . 

d4)(x)  _  ^ 
ax 


(3.2-1  la) 


(3.2-1 1b) 
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Evaluating  these  conditions  using  the  phase  given  in  Eq.  (3.2-10),  one 
critical  point.  xi,  is  found  to  be  located  at 


x  =  x  = 


(3.2-12) 


Using  a  critical  point  of  the  first  kind,  the  integral  in  Eq.  (3.2-8)  can  be 
approximated  by 


G 


•  nr^6^/(2Z„)  /  2n  \* 


ilZjpl  ♦(Xj) 


(3.2-13) 


where  the  subscript  (s)  is  used  to  denote  the  stationary  phase 
approximation  to  G(fx).  By  direct  cakulation,  the  stationary  phase 
approximation  to  the  Fresnel  zone  transform  is  given  by 


GX'>  = 


/(2Zp) 


(3.2-14) 


for  large  values  of  |Zf|. 

In  the  image  reconstruction  experiments  discussed  in  Chapt.  2  it  was 
shown  that  after  only  a  few  iterations,  the  magnitude-only 
reconstructions  improve  further  away  from  the  Fourier  transform  plane. 
From  Eq.  (3.2-14)  we  see  that  for  large  values  of  |2f1,  or  equivalently  an 
increase  in  the  distance  away  from  the  Fourier  plane,  the  magnitude  of 
the  Fresnel  zone  transform  is  approximated  by  a  scaled  version  of  the 
unknown  image.  As  the  offset  from  the  Fourier  plane  increases,  the  given 
Fresnel  zone  transform  magnitude  data  more  closely  resemble  the 
unknown  image  data. 
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3.2.3  Analytic  image  Reconstruction 

The  iterative  method  of  image  reconstruction  from  magnitude-only 
data  starts  with  an  initial  random  guess,  g^°*(x),  for  the  object  which  fulfills 
any  known  object  constraints.  This  random  guess  is  denoted  by  r(x),  the 
Fresnel  zone  transform  is  denoted  by  R(fx),  and  can  be  written  in  polar 
form  as  |R(fx)|exp{i6(fx)}.  The  algorithm  continues  by  combining  the 
known  transform  magnitude  information,  |G(fx)|,  with  a  guess  for  the 
phase  of  the  Fresnel  zone  transform.  The  starting  guess,  G^^^fx),  for  the 
Fresnel  zone  transform  is  given  by 


G‘®(4)=|G(4)|e 


ie<°V ) 

X 


(3.2-15) 


where  exp{i0<°>(fx)}  is  the  initial  guess  for  the  phase,  which  is  set  equal  to 
the  phase  of  |R(fx)|exp{i0{fx)}. 

The  first  iterative  guess  for  the  image,  g'’>(x),  is  now  obtained  by  first 
calculating  the  inverse  Fresnel  zone  transform  of  G^°>(fx)» 


= 


i  ni^2  Zjj,/ 


f  +0D 

1.  ' 


\Otft\e 


i  e‘°V ) 

X 


i  2n  f  X 


(3.2-16) 


and  then  imposing  any  known  image  constraints.  An  illustrative  example 
showing  the  result  of  this  first  iteration  is  discussed  in  Sect.  3.6. 

The  iteration  continues  by  using  g^^^x)  as  the  new  guess  for  the  image 
and  calculating  the  Fresnel  zone  transform  of  this  new  guess.  The  phase 
from  this  Fresnel  zone  transform,  exp{i0<’Hfx)},  is  then  combined  with  the 
known  magnitude,  |G(fx)|,  to  form  a  new  guess,  G<’>(fx),  for  the  Fresnel 
zone  transform.  Imposing  any  image  constraints  on  the  inverse  Fresnel 
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zone  transform  of  G<’>(fx)  yields  a  new  guess  for  the  image,  g^^*(x).  The 
algorithm  continues  to  transform  back  and  forth  while  applying  any 
known  image  constraints  and/or  any  known  transform  information. 


3.2.4  Approximation  of  Reconstructed  Image 

Using  the  approximation  to  the  Fresnel  zone  transform  developed  in 
Sect.  3.2.2  it  is  possible  to  estimate  the  quantities  needed  for  the  first 
iterative  guess  defined  by  Eq.  (3.2-16).  By  using  these  estimates,  an 
analytic  approximation  to  the  first  iterative  guess  for  the  unknown  image 
will  be  obtained.  Let  gs*’’(x)  denote  an  estimate  of  the  first  iterative  guess 
for  the  image,  g<’'(x).  This  estimate  is  then  given  by 


inz^2  Zp! 


-f  00 


ii\ft 


i2nfx 


(3.2-17) 


in  which  the  stationary  phase  approximation  to  the  known  magnitude 
information  and  the  initial  guess  for  the  unknown  phase  are  denoted  by 
|Gs(fx)|  and  exp(i0s^°>(fx)}  respectively.  Using  Eq.  (3.2-14)  as  an  estimate  for 
the  known  magnitude  information,  (G(fx)I  is  given  by 


|G  (/,)!  = 


(3.2-18) 


An  estimate  for  the  initial  guess  for  the  phase  of  the  Fresnel  zone 
transform,  exp{i6<°’(fx)},  is  given  by 


e  *  *  =  e- 


I  n  /4 


I  n  6^/(2Zp) 


(3.2-19) 
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Using  the  estimates  given  in  Eqs.  (3.2-18)  and  (3.2-19),  Eq.  (3.2-17)  can  be 
written  as 


(1) 


(x)  =  c 


2  2 
inx  2 Z^l  b 


X 


^'4 


i2nf  X 

e  *  df 


(3.2-20) 


If  this  inverse  Fresnel  zone  transform  is  viewed  as  the  inverse  Fourier 
transform  of  a  product,  which  is  equal  to  the  convolution  of  the  inverse 
Fourier  transforms,  Eq.  (3.2-20)  can  be  expressed  as 


8 


111 

$ 


(x)  = 


-ib^  .ni^2  ZJ  ^ 


X 


2Z  -,ni^2Z  /6^ 


2Z^  2Z^x 


(3.2-21) 


where  G/:(f )  is  the  Fourier  transform  of  g(x)  and  0  denotes  a  convolution. 
This  convolution  can  be  written  out  explicitly  as 


2Z^ 

V  2Z^  ) 

r  +« 

/  —  OR 

3/2  i  n  2  Zp  /  6^ 

e  ^  X 


-iT,(x-x-)^2zib^  (3.2-22) 

dx  , 
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which  when  expressed  in  the  form  of  an  inverse  Fourier  transform 
becomes 


'+«  2Z  z’  -inx'^2Z^/b^  i2n(x2Z„/6^x’ 

G/—^)  e  ^  e  ^  dx' .  (3.2-23) 

—  OB  ^ 


Again  this  inverse  Fourier  transform  of  a  product  can  be  written  as  the 
convolution  of  the  inverse  Fourier  transforms,  or 


2  2 
ix'^2  Z^i  }> 


(3.2-24) 


Writing  this  convolution  out  explicitly  as 


iZ^  2n<x-*')^ 


dx' 


(3.2-25) 


presents  the  integral  in  a  form  which  can  be  approximated  by  the  method 
of  stationary  phase.  In  this  case  the  argument  of  the  complex  exponential 
is  given  by  i|ZFl<I>(x'),  where 

4>  (x' )  =  s^nfZ^  (2  n /6^)^  z  —  x'^  (3.2-26) 

The  stationary  phase  point  is  found  by  evaluating  the  expression 

(3.2-27) 

dx' 

The  stationary  phase  point  is  located  at 
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(3.2-28) 


X 


X. 


Using  the  stationary  phase  approximation  for  the  integral  in  Eq.  (3.2-25) 
leads  to  the  approximation  for  the  first  iterative  guess  for  the  image; 

^2 

(3.2-29) 

Recall  that  gs^'^(x)  is  an  analytic  estimate  for  the  first  iterative  guess,  g^'Hx). 
By  the  approximation  given  in  Eq.  (3.2-29),  it  has  been  shown  that  for 
large  values  of  |Zf|,  the  reconstruction  after  the  first  iteration  is 
approximately  equal  to  the  unknown  image. 
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3.3  Summary  of  Results  for  Two  Dimensions 

The  results  of  the  approximations  if  using  a  two-dimensional  image  are 
directly  analogous  to  the  one-dimensional  case  discussed  above. 

Consider  a  general  function,  g(x,y),  and  its  associated  Fresnel-zone 
transform  G(fx,  fy).  The  Fresnel-zone  transform  of  g(x,y)  is  given  by 


r  +00  .  +09 

j(.xcf)K(x,y,  f^,f^)dxdy 


(3.3-1) 


By  direct  calculation,  one  can  readily  verify  the  following  inversion 


relationship 


#  4*09  f  +a» 

g(x,y)=  I 


The  Fresnel-zone  transform  kernel  K  (x,y;  fx.fy)  is  given  as 


2  2 

-ina(z  )  -i2n(f  x  +  f  y) 

Kix,y,f^.f^)^e  e  *  ^ 


(3.3-2) 


(3.3-3) 


The  conjugate  complex  is  denoted  by  *,  and  a  is  a  constant  which  is  a 
measure  of  the  distance  away  from  the  Fourier  plane  and  is  defined  in  Eq. 
(2.2-5). 

Consider  a  function  g(x,  y)  which  is  non-zero  only  over  a  finite  region, 


g{x,y)  =  Q  for  |i|  )  6  or  l^l  )  6 


(3.3-4) 
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The  number  of  Fresnel  zones,  Zp.  is  defined  as  the  number  of  multiples  of 
2n  that  the  argument  of  the  complex  Gaussian  goes  through  while 
moving  along  either  the  x  or  y  axis  from  x  =  y  =  0  to  either  x  s  b  or  y  s  b. 
The  two-dimensional  Fresnel  zone  transform  given  in  Eqs.  (3.3-1)  and 


(3.3-3)  can  be  rewritten  as 

G{fJ^)  =  e  *  >  ^ 


X 


2  2 
^  f.s 


,+br+b  ^  IV  2Z^/\  22jp/  I 

I  g(x,y)e  dxdy  (3.3-5) 

where  the  non-zero  region  of  g(x,  y)  is  reflected  in  the  bounds  of 
integration. 

The  argument  of  the  complex  Gaussian  in  this  two  dimensional 
integral  is  equal  to  i|ZF|4>(x,  y),  where  4>(x,  y)  is  given  by 

b^f  ^  b^f  ^ 

4>(xo')=  -sgntZjpI(2n/6^)[^i+ ‘  (3-3-6) 


The  stationary  phase  points  are  found  at  locations  where 

Ul  <  *  .  bl  <  b  . 


(3.3-7a) 


a<|)(x)  a<i)(>) 


(3.3-7b) 


One  stationary  phase  point  is  found  to  be  located  at 


6V  6V 

_  —  *  _  _  / 

■  2Z^’  ^"^1“  ■  2Z^ 


(3.3-8) 
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Using  the  method  of  stationary  phase  to  evaluate  the  double  integral 
in  Eq.  (3.3-5)  leads  to  the  approximation 


G 


i  n  ( 


f[  +  ^  )6^(2Zp) 


g  ( 


-b^f  -b^f 

*  X 


’  2Z^ 


(3.3-9) 


for  large  values  of  |Zf|.  This  expression  is  the  two  dimensional  form  of  Eq. 
(3.2-14). 

If  |Gs(fx,  fy)|  is  used  for  an  estimate  of  lG(fx,  fy)| ,  and  exp{i0s‘’>(fx,  fy)} 
for  an  estimate  of  exp{i6^^’(fx,  fy)}  in  a  single  iteration  of  the 
reconstruction  algorithm,  the  same  derivation  that  was  discussed  in  Sect. 
3.2.4  leads  to  the  relation 


(^)  gi^.y)  ■  (3.3-10) 

This  expression  is  directly  analogous  to  Eq.  (3.2-29)  in  Sect.  3.2.4.  For  large 
values  of  |Zf|,  the  first  iteration  of  the  reconstruction  algorithm  produces 
an  image  which  is  a  close  estimate  to  the  unknown  image. 


( 


4 


-58- 


3.4  Geometrical  Interpretation  of  Stationary  Phase 

it  is  possible  to  derive  the  location  of  the  stationary  phase  point  given 
in  Eq.  (3.2-12)  using  an  argument  based  upon  the  principles  of 
geometrical  optics.  The  expression  for  the  number  of  Fresnel  zones,  Zf, 
can  be  written 


fe2(i_d/F)i  (3.4-1) 

z  = _ - _ 

^  1  F  X  ■ 

In  the  geometrical  optics  limit,  i.e.  as  the  wavelength  x  goes  to  zero,  the 
number  of  Fresnel  zones  goes  to  infinity  (assuming  di^F).  So  the 
requirement  for  the  stationary  phase  approximation,  that  the  number  of 
Fresnel  zones  be  very  large,  is  equivalent  to  having  the  wavelength  go  to 
zero. 

■^he  method  for  performing  an  optical  Fresnel  zone  transform  is  shown 
in  Fig.  3-1  with  a  geometrical  ray  originating  at  a  point  Xo  and 
propagating  parallel  to  the  optical  axis.  After  being  refracted  by  the  lens 
of  "'ocal  length  F,  this  ray  intersects  the  u,v-plane  (i.e.  the  Fresnel  zone 
trai  sform  plane)  at  a  height  Uq.  By  similar  triangles,  the  relationship 

fo  ^  (3.4-2) 

F  F-dj 

is  obtained.  This  relationship  can  be  solved  for  x© ,  to  give 

(3.4-3) 

*0  1  -  d^/F 

Using  the  definition  for  the  number  of  Fresnel  zones  given  by  Eq.  (3.4-1) 
and  the  generalized  frequency  variable, 


U 


an  expression  for  the  stationary  phase  point,  x©,  is  given  by 


(3.4-4) 
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(3.4-5) 


f  \F 

'  X 

X  =  —  , 

^  2Z^\F/h^ 

or 

_fy  (3.4-6) 

2Z^ 

This  is  equivalent  to  the  stationary  phase  point  given  by  Eq.  (3.2-12). 


"Geometrical  optics  ray" 


Fig.  3-1  Optical  Fresnel  zone  transform  system  showing  geometrical  optics  ray. 


In  the  asymptotic  expansion  of  the  Fresnel  zone  transform,  the  leading 
term  is  attributed  to  stationary  phase  points.  The  location  of  these 
stationary  phase  points  are  the  same  as  the  points  obtained  by  tracing  an 
undiffracted  geometrical  optics  ray  through  the  imaging  system.  The 
existence  of  a  stationary  phase  point  is  easily  verified  by  tracing  a 
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geometrical  optics  ray  through  the  system,  rays  which  are  not  obstructed 
give  rise  to  stationary  phase  points. 
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3.5  Critical  Points  of  the  Second  Kind 

Section  3.3  discussed  critical  points  of  the  first  kind,  or  stationary  phase 
points.  The  stationary  phase  approximation  given  in  Eq.  (3.3-9)  is  only  the 
first  term  in  an  asymptotic  expansion  of  the  Fresnel  zone  transform 
integral.  This  section  will  discuss  critical  points  of  the  second  kind,  which 
constitute  a  higher  order  term  in  the  asymptotic  expansion.  Although 
there  is  a  close  similarity  between  one-  and  two-dimensional  integrals  for 
critical  points  of  the  first  kind,  there  is  not  such  an  easy  transition  for 
critical  points  of  the  second  kind.  For  this  reason  only  the  case  of  a  two- 
dimensional  image  will  be  discussed  here. 

Consider  once  again  the  function  given  in  Eq.  (3.3-6); 


2z) 

Critical  points  of  the  second  kind  are  located  along  the  boundary  curve  of 
the  region  of  integration.  If  a  displacement  along  the  boundary  curve  is 
denoted  by  dl,  a  critical  point  of  the  second  kind  is  located  where  the 
function  <l>(x,y),  is  stationary  along  the  boundary  curve,  or 

3<|>(j0')  _  ^  (3.5-2) 

dl 

A  rectangular  region  with  |x|  ^  b  and  |y|  ^  b  is  shown  in  Fig.  3-2.  The 
calculation  of  the  critical  points  is  now  done  by  finding  points  along  the 
rectangular  boundary  where 
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Fig.  3-2  Rectangular  region  of  integration  in  the  x,y  plane. 


3/  ”  3x  dl  dy  dl~ 


(3.5-3) 


Considering  the  direction  of  travel  to  be  counterclockwise  along  the 
boundary  of  integration,  the  partial  derivatives  are  calculated 


3  X  dy 

up  the  right  hand  side:  —  =  0  ,  —  =  +1, 

dl  dl 


(3.5-4a) 


across  the  top:  —  =  —  1  ,  —  =  0, 

dl  dl 


(3.5-4b) 


3  I  3^ 

down  the  left  hand  side:  —  =  0  ,  —  =  —  1, 

dl  dl 


(3.5-4C) 
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and  across  the  bottom: 


(3.5-4d) 


d  X  ^ 

=  +1  ,  =  0 
31  31 


The  remaining  partial  derivatives  are  given  by 


3<|)(x^y)  2n  ®  'x 


(3.5-5a) 


34*  (*i^)  2n  ®  'y 


(3.5-5b) 


Each  side  is  found  to  contain  one  critical  point  of  the  second  kind.  These 
four  points  are  located  at; 


X  =  +  6  ,  V  s=  —  ■  . 

a  '  2Z/ 


(3.5-6a) 


*  2Z^  * 


(3.5-6b) 


_  _ ,  _  _ _ 

’  ^c~  2Z^  ’ 


(3.5-6C) 


2Zp 


(3.5-6d) 


The  Fresnel  zone  transform  given  in  Eq.  (3.3-5)  is  now  approximated  by 
the  first  two  terms  in  an  asymptotic  series  as 


(3.5-7) 
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where  Gs(fx.  fy)  is  due  to  critical  poir^ts  of  the  first  kind  and  is  given  in  Eq. 
(3.3-9).  G2(fx.  fy)  is  equal  to  the  sum  of  the  contributions  from  each  of  the 
four  critical  points  of  the  second  kind.  The  contribution  from  each  of 


exp{+in/4}  fi>0  (3.5-10) 

exp{~  in/4}  n<  0 


and  R  is  the  radius  of  curvature  (which  is  infinite  for  this  rectangular 
region  of  integration). 

Evaluating  Eq.  3.5-8  at  each  critical  point  and  summing,  gives 
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where 


Q  = 


(2n)^i 

|Z^|^(2in^ 


sgviZj^  exp{  +i  n  /  4  }  . 


(3.5-12) 


(3.5-13) 


Note  that  these  terms,  which  are  due  to  the  contributions  of  the  critical 
points  of  the  second  kind,  are  of  the  order  The  term  due  to  critical 

points  of  the  first  kind  goes  as  (Zf)  ’.  For  a  very  large  number  of  Fresnel 
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zones,  the  Fresnel  zone  transform  is  dominated  by  the  term  Gs(fx,  fy), 
which  is  due  to  critical  points  of  the  first  kind. 

In  Eq.  3.3-9  it  is  very  easy  to  see  the  behavior  of  the  magnitude  and 
phase  of  the  Fresnel  zone  transform  due  to  the  critical  points  of  the  first 
kind.  It  is  not  so  easy  to  understand  the  structure  of  the  Fresnel  zone 
transform  due  to  critical  points  of  the  second  kind.  One  feature  that  can 
be  seen  readily  is  the  impact  of  the  rectangular  shape  of  the  image  on  the 
transform.  By  looking  along  the  fx  axis,  i.e.  by  letting  fy  =  0,  the 
contributions  from  the  points  (xa,  ya)  and  (xc,  yc)  depend  only  on  the 
image  points  at  x=  ±b  and  y  =  0,  the  contributions  due  to  the  critical 
points  (xb,  yb)  and  (x^,  yd),  are  proportional  to  the  magnitude  of  the 
image  along  the  edges  y=  +  b  and  y  =  -b.  Thus,  critical  points  of  the 
second  kind,  which  lie  at  the  edge  of  the  image,  yield  contributions  to  the 
Fresnel  zone  transform  which  contain  information  about  the  edges  of  the 
image. 
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3.6  Illustrative  Example 

As  described  in  Sect.  2.2,  the  numerical  computation  of  a  Fresnel  zone 
transform  is  done  using  a  fast  Fourier  transform  (FFT)  routine.  The  object 
to  be  transformed  is  multiplied  by  the  complex  Gaussian  corresponding  to 
a  particular  Fresnel  zone;  the  FFT  of  this  product  is  equal  to  the  Fresnel 
zone  transform.  Although  the  stationary  phase  approximation  discussed 
above  improves  as  the  parameter  |Zf|  increases,  the  sampling 
requirements  of  the  numerical  calculation  of  the  Fourier  transform  place 
an  upper  limit  on  the  practical  size  of  |Zf|.  A  discussion  of  sampling 
considerations  is  given  in  detail  in  Sect.  2.2.3. 

Figure  3-3  shows  the  original  digitized  image  of  a  sailboat  that  was 
used  in  this  example  of  magnitude-only  image  reconstruction.  In  this 
experiment  the  256x256  image  is  padded  with  zeros  to  give  a  512x512 
image.  For  display  purposes.  Fig.  3-3  shows  only  the  central  256x256 
pixels. 

The  phase  and  magnitude  of  a  Fresnel  zone  transform  of  the  image 
shown  in  Fig.  3-3  are  shown  in  Figs.  3-4(a)  and  (b)  respectively.  Figures  3- 
4(a)  and  (b)  are  each  displayed  with  512x512  pixels.  This  Fresnel  zone 
transform  was  calculated  for  16  Fresnel  zones.  Note  that  both  the  phase 
and  magnitude  have  structures  which  are  similar  to  that  predicted  by  Eq. 
(3.3-9).  The  phase  of  the  Fresnel  zone  transform  shown  in  Fig.  3-4(a)  is 
dominated  by  the  ring-like  structure  associated  with  a  complex  Gaussian. 
The  magnitude  of  the  Fresnel  zone  transform  shown  in  Fig.  3-4(b)  has  a 
rough  outline  of  the  sails  from  the  sailboat,  and  the  light  sky  and  dark  sea 
are  transposed.  This  general  structure  of  the  magnitude  of  the  Fresnel 
zone  transform  was  previously  predicted  and  can  be  verified  by  simply 
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Fig.  3-3.  Original  digitized  image  used  in  the  reconstruction 
experiments. 

rotating  Fig.  3-4(b)  by  180°;  i.e.,  as  shown,  the  sailboat  is  upside  down 
and  reduced  in  size. 

In  the  iterative  method  of  reconstruction,  the  given  magnitude 
information  shown  in  Fig  3-4(b)  is  multiplied  by  the  phase  from  the 
Fresnel  zone  transform  of  a  random  guess  for  the  unknown  image.  An 
inverse  Fresnel  zone  transform  then  gives  the  first  iterative  guess  that  is 
shown  in  Fig.  3-5(a).  For  display  purposes,  only  the  central  256x256 
pixels,  out  of  a  field  of  512x512  pixels,  are  shown.  In  this  first  iterative 
reconstruction  the  sailboat  is  easily  recognized,  and  even  the  small  details 
of  the  lighthouse  can  be  discerned.  This  example  is  a  demonstration  of 
the  fact  that  the  reconstruction  after  the  first  iteration  is  a  very  good 
approximation  to  the  original  image,  as  was  predicted  by  the  relationship 
given  in  Eq.  (3.3-10).  For  comparison  with  the  reconstructions  discussed  in 
Chapt.  2,  the  iterative  reconstruction  after  five  iterations  is  shown  in  Fig. 


(a)  Phase 


(b)  Magnitude 


Fig.  3-4.  {a)Fresnel  zone  transform  phase,  at  16  Fresnel  zones,  of  the  image  shown  in 
Fig.  3-3  The  phases  from  -  n  to  +  n  are  shown  as  grey  levels  from  black  to 
white.  (b)Fresnel  zone  transform  magnitude,  at  16  Fresnel  zones,  of  the 
image  shown  in  Fig.  3-3  The  magnitude  is  shown  on  a  linear  scale  such  that 
|G(fx,  fy)|i  =  0  is  black  and  the  maximum  value  of  (G(fx,  fy)|i-  is  white. 
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Chapter  Four 


The  Importance  of  Edges  in  Phase-only  Reconstructions 

4.1  Introduction 

The  importance  of  the  phase  information  in  image  formation  and 
signal  processing  has  been  discussed  in  the  literature  [see  for  example, 
O'Neil  and  Walther,  1963;  Huang  et  al,  1975;  Oppenheim  et  al,  1979; 
Oppenheim  and  Lim,  1981]  The  Kinoform,  a  device  which  operates  only 
on  the  phase  of  an  incident  wave,  has  been  developed  [Lesem  et  al,  1969]. 
Kermisch  [1970]  has  given  a  quantitative  analysis  of  image  reconstruction 
from  phase-only  data  for  the  case  of  a  diffusely  reflecting,  coherently 
illuminated  surface.  For  the  case  of  Fourier  transform  data  compression, 
Pearlman  and  Gray  [1978]  have  shown  that  the  real  and  imaginary  parts  of 
the  data  should  be  encoded  with  equal  precision,  but  that  in  polar  form 
the  phase  angle  is  more  important  than  the  magnitude.  Using  only  the 
phase  information  in  correlation  filters  has  also  been  investigated  [Horner 
and  Gianino,  1984;  Gianino  and  Horner,  1984] 

A  simple  reconstruction  experiment  can  be  used  to  illustrate  the 
importance  of  the  phase  information  in  the  Fourier  transform.  Figure  4- 
1(a)  is  the  original  256  X  256  X  8-bit  image  denoted  by  g(x,  y).  The  Fourier 
transform  of  this  image  is  denoted  by  G(fx,  fy),  and  is  a  complex  quantity 
with  both  a  magnitude  and  phase.  This  complex  quantity  can  be 
represented  in  polar  form  as 

=  ,  (4.1-1) 
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where  |G(fx,  fy)|  is  the  magnitude  of  the  transform,  and  exp{i6(fx,  fy)}  is 
the  exponential  phase  of  the  transform. 

The  resulting  image  of  a  phase-only  inverse  transform  is  shown  in  Fig. 
4-2(b).  This  image  is  obtained  by  taking  the  inverse  Fourier  transform  of 
the  quantity  exp{i6(fx,  fy)}.  The  magnitude  information  has  been 
discarded  by  setting  the  magnitude  equal  to  unity  in  the  transform  plane. 
The  fact  that  the  edges  are  the  significant  features  which  are  retained  in 
this  experiment  is  easily  seen.  This  type  of  experiment  was  first  reported 
by  Oppenheim,  Lim,  Kopec,  and  Pohlig  [1979]. 

One  idea  which  would  explain  this  phenomenon  is  based  upon  the 
technique  of  edge  enhancement  through  high-pass  filtering.  For  a  wide 
class  of  interesting  images,  the  magnitude  of  the  Fourier  transform 
generally  behaves  as  1/r  for  r^O,  where  r2  =fx2  +  fy2  [Scott,  1963;  Kinziy 
etal,  1968).  Dividing  both  sides  of  Eq.  (4.1-1)  by  |G(fx,  fy)|  gives 


(4.1-2) 


■x-y 


Filtering  the  Fourier  transform  by  |G(fx,  fy)|~’  amplifies  the  high 
frequencies  and  attenuates  the  low  frequencies.  In  the  image  domain  this 
has  the  effect  of  enhancing  the  edges. 

While  setting  the  magnitude  equal  to  unity  leads  to  a  very  edge 
enhanced  image,  it  is  also  interesting  to  combine  the  phase  from  one 
image  with  the  magnitude  of  a  second  image.  In  this  case,  the  re*,  jiting 
image  contains  many  of  the  features  of  the  image  from  which  the  phase 
was  taken  [Oppenheim  and  Lim,  1981].  The  features  which  are  preserved 
are  again  the  high-contrast  edges  of  the  image. 
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(b)  Phase-only  Inverse  Transform 


Fig.  4-1 


(a)  Original  digitized  image  used  in  this  reconstruction  demonstration. 

(b)  Magnitude  of  the  phase-only  reconstruction  obtained  by  setting  the  Fourier 
transform  magnitude  equal  to  unity  and  calculating  an  inverse  Fourier 
transform. 
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Section  4.2  contains  a  discussion  of  two  examples  of  phase-only 
reconstructions  which  can  be  done  analytically.  These  examples 
demonstrate  the  preservation  of  the  edges  in  the  Fourier  transform  phase 
information. 

Reconstructions  from  Fresnel  zone  transform  phase  information  have 
been  shown  in  Sect.  2.3.  There  is  a  slight  degradation  in  the  quality  of 
these  reconstructions  as  the  number  of  Fresnel  zones  is  increased.  Section 
4.3  will  present  two  analytic  reconstructions  from  Fresnel  zone  transform 
phase-only  data.  The  approximations  used  in  these  examples  will  assume 
that  the  number  of  Fresnel  zones  is  small.  For  a  small  number  of  Fresnel 
zones,  the  phase  information  in  the  Fresnel  zone  transform  preserves 
information  about  the  edges  in  the  image. 
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4.2  Analytic  Fourier  Transform  Phase-only  Reconstructions 

The  rectangle  function,  shown  in  Fig.  4-2,  is  a  simple,  non-negative 


2 

4>  ^ 

o 
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Fig.  4-2  Rectangle  function  denoted  by  rectI(x-xo)/bl. 

function  of  finite  extent  which  has  sharp  edges.  This  function  is  denoted 
by: 


(4.2-1) 


The  Fourier  transform  of  g(x),  G(f),  is  given  by 

_  «fi(n  6  /) 

Gif)  =  exp{-i2nfx^b  . 

This  can  be  written  as  a  magnitude. 


icv,l=  1 

n  I/"! 


and  a  phase 


1 

f-T: 


(4.2-2) 


(4.2-3) 


(4.2-4) 


exp{iQif)}=exp{-i2nfx^sgnifl  - 1 +2rect(b/)  ^  -combf - —)  , 

2  \  2/ b  } 


where  0  denotes  a  convolution.  This  way  of  expressing  the  phase 
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function  will  prove  useful  in  the  computation  of  the  phase-only 
reconstruction. 


An  expression  for  the  phase-only  reconstruction  is  found  by  taking  the 
inverse  Fourier  transform  of  Eq.  4-6,  and  is  given  by 

The  magnitude  of  this  expression  is  shown  in  Fig  4-3.  Note  that  the 


Fig.  4-3  Magnitude  of  the  Fourier  transform  phase-only  reconstruction  (solid  line)  and 
the  original  rectangle  function  (dashed  line). 

Structure  is  dominated  by  the  presence  of  large  spikes  at  x  =  xq  ±  b/2,  the 
location  of  the  edges  in  the  original  rectangle  function.  There  is  no  spike 
at  the  center  of  the  sinc{(x-xo)/b}  envelope,  and  the  higher  order  spikes 
are  attenuated  by  this  envelope  function. 

As  a  second  example  which  shows  the  importance  of  edges  in  the 
phase  data,  consider  the  hat  function  shown  in  Fig  4-4.  This  function  is 
generated  by  adding  two  rectangle  functions  together,  as 


IH 
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Hat  function 
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Fig.  4-4  Sum  of  rectangle  functions,  producing  the  hat  function. 


(4.2-6) 


The  Fourier  transform  of  Eq.  (4.2-6)  can  be  written  as  a  magnitude 


(Xj-Xj) 


_ _ _  2  |«m(n6/)  „ 

|G(/^)|= - — -  CDS  2nf - - 

n  /I  12 


(4.2-7) 


and  a  phase 


1 


f  ^*1  *2^  f  b  f  2b  \ 

exp{iQ{f)}=  exp  —i2nf - - -  sgnif)  —  1 2rec<(6/)  0  -combi  - —  I  X 

I  2J  I  2  \  21b  / 


The  inverse  Fourier  transform  of  Eq.  (4.2-8)  is  given  by: 


(4.2-8) 


/Xj+Xj-,  _j  x2  ( 2x\ 

8|^z- — ^ — j  ®  •; —  ®  -8(x)  +  siruK")  -oom6^ jexp{inx/6} 


0  —  8(x)  +  sincf - ^  -  comb 

\  X,  -  Xj  /  (Xj  -  X^)  \  X„  -  X,  /  J 


2  ■*! 
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and  the  magnitude  is  shown  in  Fig.  4-5.  Notice  that  once  again,  the 
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Fig.  4-5  Magnitude  of  the  Fourier  transform  phase-only  reconstruction  (solid  line)  and 
the  original  hat  function  (dashed  line). 

general  structure  of  the  phase-only  reconstruction  is  dominated  by  the 
appearance  of  large  spikes  at  the  same  locations  of  the  edges  of  the 
original  function. 

Both  of  these  one-dimensional  examples  have  shown  how  the 
information  about  the  edges  in  an  object  is  preserved  in  the  phase  of  the 
Fourier  transform.  In  performing  a  phase-only  reconstruction,  the  edge 
itself  is  not  reconstructed,  but  rather  a  strong  discontinuity  arises  at  the 
location  of  the  edge. 
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4.3  Analytic  Fresnel-zone  Transform  Phase-only  Reconstructions 

Now  consider  the  case  of  the  rectangle  function  shown  in  Fig.  4-2 
under  the  process  of  a  Fresnel-zone  transform  phase-only  reconstruction. 
The  Fresnel-zone  transform  of  Eq.  4-3  is  written  as 


G(/)  =  l  red^ — exp(—inax^)  exp(—i2nfx)dx.  ^ 

Using  the  convolution  property  of  Fourier  transforms,  this  integral  can  be 
treated  as  the  Fourier  transform  of  a  product.  The  result  is  a  convolution 
of  Fourier  transforms  given  by. 


G(/‘)=  exp{-i2nfxQ^bsinc(bft 


®  - -  expdn  /*/a)  , 

(io)* 


which  when  written  in  integral  form  becomes 


(4.3-2) 


G(r)  = 


exp(-inax^) 

n(ia)* 


exp{-i2nrx^j  ^ exp{i -l/'-(ax^-hPl^}dr  . 

(4.3-3) 


A  useful  approximation  to  this  integral  is  possible  with  the  assumption 
that  a  is  small,  i.e.  the  Fresnel  zone  is  located  close  to  the  Fourier 
transform  plane.  It  is  then  possible  to  use  a  stationary  phase 
approximation  (see  also  Sect.  3.2)  to  evaluate  the  integral  in  Eq.  (4.3-3). 
The  stationary  phase  approximation  leads  to  an  approximate  expression 
for  the  Fresnel-zone  transform  of  rect{(x-xo)A)) ; 


G(f)‘~exp(.-inax^)  exp{—i2nfx^b 


sini  n  i»(axjj+  f)] 


n  b{ax.+  f) 


Which  can  be  separated  into  a  magnitude; 

I  I  sin(  n  6(ax  +  / )) 

|G(r)| - 

n 


a*o+  f  I 


(4.3-4) 


(4.3-5) 
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and  a  phase; 


exp{i6{f)} -  expi-ina x^)  exp{-i2nfx^sgnlaxQ+f]  X 


I  M  (4-3-6; 

Setting  the  magnitude  equal  to  unity  and  taking  the  inverse  Fresnel- 
zone  transform  of  the  phase  leads  to  the  expression  for  the  phase-only 
reconstruction : 

etp{— i2naxjj(x-ijj)}  exp{—i2nax^} 

- - -  _ -  0 

m(x— x^  ini 


/*“  *o\  2  1 


(4.3-7) 


The  magnitude  ofthis  complex-valued  expression  is  shown  in  Pig.  4-6.  The 


Pix«l  Nuaibsr 

Fig.  4-6  Magnitude  of  the  Fresnel-zone  transform  phase-only  reconstruction  (solid 
line)  and  the  original  rectangle  function  (dashed  line). 


result  is  a  structure  which  is  very  similar  to  the  Fourier  transform  phase- 
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only  reconstruction  shown  in  Fig.  4-3.  The  dominant  feature  of  the 
structure  is  a  large  spike  at  the  locations  of  the  edges  in  the  rectangle 
function. 

As  a  final  example,  consider  once  again  the  sum  of  two  rectangle 
functions  shown  in  Fig.  4-4.  Using  the  same  argument  for  the  stationary 
phase  approximation  that  was  presented  above,  the  Fresnel  zone 
transform  of  Eq.  (4.2-6)  is  approximately  equal  to 

G(f)  ~  2b  exp(— inax^)  exp{—i2nfx^  X 


«fitn6(ax»+  /)] 


cos 


(Xj-*,) 

2n— —  (A+QXg) 


n  6(aXg+  f) 

where  X3  =  (x^  -t-  X2)/2.  The  phase  part  of  this  expression  is 


«rp{i0(/)}-erp(-inaXg)  exp{-i2nfx^sgnlf+  qx^]  X 


(4.3-8) 


—  1  +  2rectlfb)  ®  -  com6l  - 

2  \  2/6 


)1 


X 


)|  .  (4.3-9) 

V  2  V  2  V2/(x,-x,)/j 


/■+^ 
2  “ 


Taking  the  inverse  Fresnel-zone  transform  of  Eq.  (4.3-9)  gives 


2  2 

«p(— inax_)  erp{+inax  )  8(x— xj 


—  exp(—i2  nax.x) 


® 


inx 


X  2 

—  8(x)  +  «nc( ~ )  -comb 
b  b 


(t) 


exp(  — i2nax-x)  cip{inx/6}|  ® 
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Fig.  4-7  Magnitude  of  the  Fresnel  zone  transform  phase-only  reconstruction  of  the 
sum  of  rectangle  functions  (solid  line)  and  the  original  function  (dashed  line). 

examples,  the  dominant  features  in  this  reconstruction  are  the  large 

spikes  found  at  the  location  of  the  edges  in  the  original  function. 

The  phase-only  reconstructions  from  Fourier  transform  data  are  real¬ 
valued  quantities.  Setting  the  magnitude  equal  to  unity  maintains  the 
property  that  the  real  part  of  the  Fourier  transform  is  even  and  the 
imaginary  part  of  the  Fourier  transform  is  odd.  A  Fresnel  zone  transform 
does  not  have  this  property,  so  setting  the  magnitude  equal  to  unity  will 
in  general  lead  to  a  complex-valued  reconstruction.  This  investigation  has 
been  limited  to  real-valued,  non-negative  images,  so  looking  at  the 
magnitude  of  the  reconstruction  is  quite  reasonable  for  both  Fourier  and 
Fresnel  transforms. 


It  has  been  shown  for  both  Fourier  transforms  and  Fresnel-zone 


transforms  located  near  the  Fourier  transform  plane,  that  the  phase  of  the 
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transform  preserves  the  information  about  the  edges  in  the  original 
image.  As  one  moves  far  away  from  the  Fourier  transform  plane,  the 
parameter  a  (and  thus,  the  number  of  Fresnel  zones)  increases.  Section  3.2 
discusses  the  effects  of  a  large  number  of  Fresnel  zones,  and  it  is  shown 
that  for  a  real-valued,  non-negative  image,  the  phase  of  the  Fresnel-zone 
transform  behaves  as  a  complex  Gaussian,  or 

erp{i0(/)}~erp(  +  in/‘^ /q)  (4.3-11) 

In  this  case  there  is  no  information  about  the  image  contained  in  the 
phase  of  the  Fresnel-zone  transform.  This  fact  accounts  for  the 
degradation  in  the  quality  of  the  phase-only  reconstructions  (for  a  set 
number  of  iterations)  from  Fresnel-zone  transform  phase-only 
information  as  the  number  of  Fresnel  zones  increases. 
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Chapter  Five 


Effects  of  Noise  on  image  Reconstruction 

5.1 .  Introduction 

The  formulation  and  techniques  used  in  image  reconstruction  from 
partial  Fresnel  zone  information  are  based  upon  properties  of  coherent 
optical  imaging,  in  coherent  imaging,  both  the  magnitude  and  phase  of  a 
signal  can  be  altered  in  a  predetermined  way  by  the  processing  system,  or 
changed  in  an  unknown  manner  by  random  fluctuations  and  noise.  In  this 
chapter  we  will  develop  a  model  for,  and  examine  the  effects  of 
coherently  added  Gaussian  noise  on  image  reconstruction  from  partial 
Fresnel  zone  information.  Sections  5.2  and  5.3  will  investigate  the  effects 
of  signal  independent  additive  noise  such  as  thermal  noise  and  amplifier 
noise,  while  Sections  5.4  and  5.5  will  examine  the  effects  of  signal 
dependent  multiplicative  noise  such  as  shot  noise,  or  film  grain  noise. 

Consider  a  scalar  component  of  the  electromagnetic  field  (dropping 
the  explicit  time  dependence)  in  which  both  the  magnitude  and  phase  are 
known.  This  signal,  Ws,  can  be  represented  in  complex  rectangular  form  as 

or  in  polar  form  as 

=  I*",!  (5.1-2) 


where 


and 


|u>J  CO«(0J  =  , 


(5.1-3) 

(5.1 -4a) 


l“.l“"(0.)  =  ^  (5.1-4b) 

The  complex  signal,  Ws.  is  shown  graphically  in  Fig.  5-1 . 


Imaginary  Axis 


Fig.  5-1.  The  true  signal,  w,  =  Xj  +  iy„  with  the  polar  representation,  Wj  = 
|wj|exp{i0,). 


When  the  noise,  n,  is  added  to  the  signal,  Ws,  we  obtain  a  noisy  signal, 
Wn,  defined  by 
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(5.1-5) 


w  =  w  +  n, 

n  $ 

where  the  noise  signal  has  been  defined  in  a  general  form  as 

n—n  +  in  . 

*  y 

The  detailed  form  of  the  probability  density  function  for  n  will  be 
discussed  in  more  detail  in  Sect.  5.2.  Writing  the  real  and  imaginary  parts 
explicitly,  we  obtain 

w^=(x^  +  n^)  +  i(y^  +  n^).  (5  -,.7) 

where  we  can  now  make  use  of  new  random  variables  Xn  and  yn  to  give 

(5.1-8) 

The  reconstructions  are  now  carried  out  using  either  the  magnitude  or 
phase  of  Wn  as  the  known  information.  Since  n*  and  ny  are  random 
variables,  wn  is  a  random  variable,  and  both  the  magnitude  and  phase  of 
Wn  are  random  variables.  One  calculation  of  interest  will  be  to  determine 
the  probability  density  functions  for  the  magnitude  and  phase  of  Wn. 

Adding  noise  to  both  the  real  and  imaginary  parts  of  the  signal  before 
calculating  either  the  magnitude  or  phase  of  the  resultant  noisy  signal  has 
several  advantages.  By  adding  the  noise  before  calculating  the 
magnitude,  we  are  assured  that  the  resultant  noisy  data,  |wn|,  is  a  quantity 
greater  than  or  equal  to  zero.  If  noise  were  added  directly  to  the  true 
magnitude  data  one  would  have  to  use  a  noise  model  which  assures  the 
non-negativity  of  the  noisy  data.  A  Gaussian  probability  density  function 
will  not  work  in  this  case.  Since  the  phase  of  the  signal  is  calculated  using 
the  ratio  of  the  imaginary  part  of  the  signal  to  the  real  part  of  the  signal. 
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any  measurement  which  records  these  quantities  will  produce  a  noisy 
estimate  of  each  of  these  values.  By  using  this  coherently  additive  noise 
model  we  are  able  to  examine  the  effects  of  noise  on  reconstructing  an 
image  both  from  magnitude-only  data  and  from  phase-only  data,  where 
the  simulated  noisy  data  is  the  result  of  a  single  noise  model. 


-90- 


5.2.  Coherent,  Additive,  Signal-independent  Noise 
5.2.1.  Graphical  Representation 

In  image  reconstruction  experiments,  either  the  magnitude,  |ws|,  or  the 
phase,  exp{i6s},  of  the  Fresnel  zone  transform  is  used.  In  this  model,  a 
noise  signal,  n,  will  be  added  to  the  true  signal,  Ws,  before  the  magnitude 
and  phase  are  calculated.  The  noise  signal  will  be  represented  as; 


n  =  X  +  iy  , 


(5.2-1) 


where  x  and  y  are  independent,  uncorrelated,  Gaussian  distributed, 
random  variables  with  zero  mean  and  variance  o2.  The  joint  probability 
density  function  fxy(x,  y),  for  x  and  y  is  given  by 


■xy 


2na  ‘ 


«P 


2a  ^ 


(5.2-2) 


where  is  the  variance  of  the  noise,  and  o  is  the  standard  deviation  of 
the  noise. 

One  realization  of  the  noise  model  discussed  above  is  shown  in  Fig.  5-2. 
The  resultant  signal,  Wn,  is  obtained  by  adding  a  noise  vector,  n,  to  the 
true  signal,  Ws.  Different  realizations  of  the  noise  result  in  the  final  signal 
lying  in  a  "noise  cloud"  about  the  true  signal.  For  the  case  shown  in  Fig.  5- 
2  the  standard  deviation  of  the  noise  is  much  less  than  the  magnitude  of 
the  true  signal,  or  o<<|ws|.  This  is  the  case  for  a  large  signal  to  noise 
ratio.  As  expected  for  the  case  of  a  large  signal  to  noise  ratio,  the 
measured  magnitude  and  phase  will  most  likely  be  very  close  to  the  true 
unknown  data.  The  variance  of  the  measured  signal  will  be  small 
compared  to  the  magnitude  of  the  signal. 
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Fig.  5-2.  Noisy  signal,  Wn,  results  from  adding  noise  vector,  n,  to  the  true  signal,  w,. 

For  the  case  shown  here  the  standard  deviation  of  the  noise  is  much  less  than 
the  magnitude  of  the  true  signal,  resulting  in  a  small  'noise  cloud’. 


The  case  of  a  small  signal  to  noise  ratio  is  shown  in  Tig.  5-3.  In  this 
instance,  o>>|w$|,  and  the  measured  data  may  differ  significantly  from 
the  true  signal.  The  variance  of  the  measured  signal  will  depend  strongly 
on  the  variance  of  the  noise,  and  not  on  the  magnitude  of  the  true  signal. 
The  true  signal  is  buried  in  the  noise.  For  the  special  case  when  |ws|  =  0., 
the  probability  density  function  for  the  magnitude  and  phase  are  the  well 
known  Rayleigh  and  uniform  distributions  respectively.  These 
distributions  will  be  shown  to  be  the  limiting  cases  as  |ws|  tends  to  zero. 
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Imaginary  Axis 


Fig.  5-3.  Noisy  signal,  Wn,  results  from  adding  noise  vector,  n,  to  the  true  signal,  Wj. 

For  the  case  shown  here  the  standard  deviation  of  the  noise  is  much  greater 
than  the  magnitude  of  the  true  signal,  resulting  in  a  large  'noise  cloud*. 


This  signal  incJependent  noise  model  has  the  property  that  the  variance 
of  the  measured  data  will  be  inversely  proportional  to  the  ratio  |ws|:o. 
Although  both  the  magnitude  and  phase  will  depend  on  the  ratio  |ws|:o, 
the  resultant  value  of  |wnl  will  be  independent  of  the  phase  angle  of  Ws, 
and^'  e  measured  phase  will  depend  upon  the  phase  angle  of  Ws  only  as  a 
bias  term.  It  will  also  be  shown  that  the  probability  density  function  for 
the  measured  magnitude  is  non-zero  only  for  magnitudes  greater  than  or 


equal  to  zero,  and  that  the  probability  density  function  for  the  phase 
angle  is  symmetric  about  the  true  phase  angle. 

5.2.2.  Analytic  Calculation 

There  is  a  certain  type  of  rotational  invariance  in  the  distributions  of 
|wn|  and  6n.  Although  the  resultant  value  for  wn  depends  strongly  on  the 
relative  sizes  of  |ws|  and  o,  the  density  function  for  |vvn|  is  independent  of 
6s,  and  the  density  function  for  6n  is  centered  about  6$.  To  simplify  the 
analytic  calculation  of  the  desired  density  functions,  we  will  consider  a 
signal  for  which  6$  is  zero.  Let 

(5.2-3) 

w  =  s.  '  ' 

« 

Consider  the  probability  density  function  fxy(x,  y)  defined  in  Eq.  (5.2-2). 
This  density  function  represents  the  "noise  cloud"  shown  in  Figs.  5-2  and 
5-3.  For  the  given  signal,  Ws  =  s,  this  noise  cloud  will  be  centered  about  the 
point  (xn  =  s,  yn*0).  In  effect  this  probability  density  function  for  the 
noise  undergoes  a  coordinate  transformation  such  that 

x  =  6oos6  — 8  ,  (5.2-4a) 

and 

y  =  bsinQ  .  (5.2-4b) 

where 

b  =  \wj  (5.2-5) 

is  the  resultant  magnitude,  and 

0  =  0^  (5.2-6) 
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is  the  resultant  phase  angle.  [For  a  similar  discussion  of  this  type  of 
calculation  see;  Papoulis,  (1965)  page  195  Ex.7<8a,  page  229  Prob.  7-4,  and 
the  discussion  leading  to  Eq.  14-63.  also;  Goodman,  (1985)  §  2.9.4] 

This  coordinate  transformation  on  the  probability  density  function 
leads  to  the  new  probability  density  function,  fbe(b,  6),  which  is  the  joint 
probability  density  function  for  the  magnitude,  b,  and  the  phase  angle,  6, 
and  is  given  by 

where  |  J  |  is  the  Jacobian  of  the  transformation  and  is  defined  as 

a*/a6  ax/ae  (5.2-8) 

dy/86  dy/ae 

Using  Eq.  (5.2-4)  to  evaluate  this  Jacobian  we  obtain 

oosO  —bsind  (5.2-9) 

=  b  . 

sinO  bcosO 

Substituting  Eqs.  (5.2-4)  and  (5.2-9)  into  Eq.  (5.2-7),  an  expression  for  the 
joint  probability  density  function  is  given  by: 

^  +  ,  (5.2-10) 

2na^  I  J 

6^0  and,  —  n  <  0  s  n  . 
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5.2.3.  Probability  Density  Function  for  Magnitude  Information 

To  calculate  the  density  function  for  the  magnitude  b,  denoted  by 
fb(b),  we  must  evaluate  the  integral  of  fb6(b.  6)  from  zero  to  2n: 


r2n 

=  I 


(5.2-11) 


Using  Eq.  (5.2-10)  this  integral  can  be  expressed  as 


(6)  =  — -  exp  I  —  (5*  +  s?)  /  2o^|  —  j  exp  sbco^  ^  ^ 

The  probability  density  function  for  the  magnitude  b  is  given  by  [see; 
Abromowitz  and  Stegun,  (1972)  Eq.  9.6.16] 


exp  -(A^  +  s^) /2o2}  6^0, 


(5.2-13) 


where  loO  is  the  modified  Bessel  function  of  the  first  kind,  defined  by  the 
power  series  [see;  Abromowitz  and  Stegun,  (1972)  Eq.  9.6.12] 


v*)=l 


(5.2-14) 


„t'o  22"(a!)=‘  ■ 


In  the  limit  as  the  true  signal,  s,  goes  to  zero,  the  probability  density 
function  for  the  magnitude  reduces  to 


(6)  =  -^  «p  -  /  2o^ 

O  ^ 


(5.2-15) 


which  is  the  well  known  Rayleigh  density  function. 

In  the  limiting  case  of  a  very  large  signal  to  noise  ratio  (i.e.  s>  >o)  the 
resultant  magnitude,  b,  will  be  approximately  equal  to  the  original  signal, 
s.  To  examine  the  behavior  of  the  probability  density  function  for  the  case 
of  a  large  signal  to  noise  ratio,  rewrite  Eq.  (5.2-13)  as 


I 


/^(6)=—  «p  I -(6- s)*  /  2o*|  ^  expj- 86*/ o*  I  ^5.2-16) 

For  s>>o  the  quantity  (sb)/a2  tends  to  infinity  for  non-zero  b.  The 
asymptotic  form  of  the  modified  Bessel  function  is  given  by 

/  (5.2-17) 

®  (2n*)* 

An  expression  for  the  large  signal  to  noise  ratio  limit  is  obtained  by  using 
this  asymptotic  expression  in  Eq.  (5.2-16),  and  is  given  by: 

f^(b)=— —  ap[ -(6-8)*  /2o* 
o  (2n)* 

Since  s>  >o,  the  exponential  term  will  have  appreciable  value  only  when 
Ssb,  and  the  term  (b/s)i  will  be  very  close  to  unity.  We  note  that  in  this 
case  the  probability  density  function  becomes  a  very  narrow,  almost  delta¬ 
like,  Gaussian  distribution  with  a  mean  equal  to  the  true  value,  s,  and  a 
variance  equal  to  the  variance  of  the  noise. 

Figure  5-4  shows  a  plot  of  ofb(b)  vs.  b/o  for  various  values  of  the 
parameter  k  =  s/o.  These  plots  show  the  smooth  transition  from  a  Rayleigh 
distribution  when  k  s  0  to  a  Gaussian  type  distribution  as  k  increases.  The 
relatively  constant  width  and  height  of  the  distribution  functions  for 
various  value  of  k  is  the  result  of  scaling  both  the  horizontal  and  vertical 
axes  by  the  standard  deviation  of  the  noise,  o. 

Three  example  histograms  and  their  corresponding  theoretical 
probability  density  functions  are  shown  in  Fig.  5-5.  For  each  each  of  these 
examples  the  known  signal  has  a  value  of  10.0  (i.e.  s=10.0),  and  the 


mr 


(5.2-18) 
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histograms  of  10000  data  points  with  noise  variance  (o2)  equal  to  1 .0, 25.0. 
and  64.0,  are  shown  in  Figs.  5-5(a),  (b).  and  (c),  respectively.  Each  of  the 
histograms  represent  the  result  of  adding  signal  independent,  Gaussian 
vvhite  noise  with  the  specified  variance  and  zero  mean  to  10000  pixels  of 
constant  value,  and  then  calculating  the  magnitude  of  the  resultant 
pixels.  Recall  that  in  the  Gaussian  noise  model  developed  here,  the  noise 
has  both  a  real  and  imaginary  parts  which  are  uncorrelated.  The  plots 


Fig.  5-4.  Probability  density  function  of  the  measured  magnitude  b  resulting  from 
noise  with  variance  added  to  a  true  signal  of  magnitude  s.  The  parameter 
k  =  s/o  is  a  measure  of  the  signal  to  noise  ratio. 
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Fig.  5-5  Predicted  probability  density  functions  and  actual  histograms  of  the  measured 
magnitude,  b,  for  (a)  s  =  1 0.0,  o  =  1 .0  (b)  s  =  1 0.0,  o  =  5.0,  and  (c)  s  =  1 0.0,  o  =  8.0 . 
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shown  in  Fig.  5-5  indicate  excellent  agreement  between  the  theoretically 
predicted  probability  density  functions  and  the  actual  histograms.  It 
should  also  be  noted  that  the  vertical  scale  is  different  for  all  three  graphs, 
and  that  the  horizontal  scale  in  Fig.  5-5(a)  has  only  one-half  the  range  of 
the  horizontal  scales  used  in  Figs.  5-5(b)  and  (c).  These  different  scales 
were  used  so  that  a  comparison  between  the  theoretical  probability 
density  functions  and  the  histograms  could  be  done  visually. 

As  the  ratio  of  the  true  signal,  s.  to  the  standard  deviation  of  the  noise, 
o,  increases  from  zero.the  probability  density  function,  fb(b),  for  the 
magnitude,  b,  changes  from  a  wide,  low  Rayleigh  distribution  to  a  very 
narrow  Gaussian. 


5.2.4.  Probability  Density  Function  for  Phase  Angle 

The  density  function  fe(6)  is  obtained  by  evaluating  the  integral 


e)  db. 


(5.2-19) 


Using  Eq.  (5.2-10)  this  integral  can  be  expressed  as 


/■g(0)=  — ^expj-s*/2o*j  I  6exp  j- 6*/2o*  + 86co«e/o*j  d6.  (5  2-20) 

The  probability  density  function  for  the  phase  angle  0  Is  given  by  [see; 
Gradshteyn  and  Ryzhik,  (1980)  Eq.  3.462.5] 


eep{— s^/2o^)  scosB  f  ^  2  11-  -  scosG 


^  =  9 

2  n 


2o(2n)* 


exp  I  —  (85in0)  / 2o  [1+2  erf\ 


II- 


\  o 


)l 


—  n  <  0  S  n  , 


(5.2-21) 
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where  erfO  is  the  error  function  defined  by 


erfix) 


=— r 

(2n)*  Jo 


exp(—u  */  2)  du  . 


(5.2-22) 


In  the  limit  as  the  true  signal,  s,  goes  to  zero,  the  probability  density 
function  for  the  phase  reduces  to 


— n  <  0  S  n  , 


(5.2-23) 


which  is  the  well  known  uniform  density  function. 

In  the  limiting  case  of  a  large  signal  to  noise  ratio  (i.e.  s>>o),  the 
resultant  phase  angle,  6,  will  be  approximately  equal  to  the  original  phase 
angle,  6s  s  0.0.  As  the  quantity  s/o  goes  to  infinity  we  obtain  an  expression 
for  the  large  signal  to  noise  ratio  limit: 

8  CD50  I  9  9 

L  (0)  =  - 7  «P  I  -  (s«m0)^  /2o  *  .  . 

®  o(2n)*  I  (5.2-24) 

The  exponential  will  be  non-zero  only  near  sin(6)x0.  (i.e.  0x0.).  Near 
8x0.,  small  angle  approximations  cos(6)  x  1.0  and  sin(6)  x6  can  be  used. 
Equation  (5.2-24)  can  now  be  expressed  in  terms  of  a  delta  function  by 
using  the  relation 


m=  „  Nexpi-N^nt^)  . 


(5.2-25) 


to  give 


Ag(0)=  6(0) 


—  n  <  0  s  n  , 


(5.2-26) 


So  for  the  case  of  a  large  signal  to  noise  ratio,  the  probability  density 
function  becomes  a  delta  function  located  at  the  true  phase  angle. 

Figure  5-6  shows  a  plot  of  fe(6)  vs.  6  for  various  values  of  the  parameter 
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Fig.  5-6  Probability  density  function  of  the  measured  phase  angle  0  resulting  from  noise 
with  a  variance  being  added  to  a  true  signal  with  a  magnitude  equal  to  s,  and 
a  phase  angle  equal  to  zero.  The  parameter  k  =  s/o  is  a  measure  of  the  signal  to 
noise  ratio. 

kss/a.  These  plots  show  the  smooth  transition  from  a  uniform 
distribution  when  k  =  0.,  to  a  Gaussian  type  distribution,  and  finally  to  a 
delta  function  as  k  tends  to  infinity. 

Three  example  histograms  and  their  corresponding  theoretical 
probability  density  functions  are  shown  in  Fig.  5-7.  For  each  each  of  these 
examples  the  known  signal  has  a  value  of  10.0  (i.e.  sslO.O),  and  the 
histograms  of  10000  data  points,  with  noise  variance  (o2)  equal  to  1.0, 
25.0,  and  64.0,  are  shown  in  Figs.  5-7(a),  (b),  and  (c),  respectively.  Each  of 
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Hi«tograB  and 
Donaity  Funotlon 


Fig.  5-7  Predicted  probability  density  functions  and  actual  histograms  of  the  measured 
phase  angle,  0,  for  (a)  s  =  1 0.0,  o  =  1 .0  (b)  s  =  1 0.0,  o  s  5.0,  and  (c)  s  s  1 0.0,  o  =  8.0 
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the  histograms  represent  the  result  of  adding  signal  independent, 
Gaussian  white  noise  with  the  specified  variance  and  zero  mean  to  10000 
pixels  of  constant  value,  and  then  calculating  the  phase  angle  of  the 
resultant  pixels.  The  plots  shown  in  Fig.  5-7  indicate  excellent  agreement 
between  the  theoretically  predicted  probability  density  functions  and  the 
actual  histograms. 

As  the  ratio  of  the  true  signal,  s,  to  the  standard  deviation  of  the  noise, 
o,  increases  from  zero,the  probability  density  function,  fe(6),  for  the  phase 
angle,  6,  changes  from  a  uniform  distribution,  to  a  narrow  Gaussian,  and 
finally  to  a  delta  function. 
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5.3  Effects  of  Additive  Noise  on  Image  Reconstruction 
5.3.1.  Magnitude-only  Reconstruction 

The  noise  model  developed  in  Sect.  5.2  was  applied  to  various  Fresnel 
zone  transforms  before  calculating  the  magnitude  data.  By  attempting  to 
reconstruct  an  image  from  this  noisy  data,  effects  of  noise  and  the  stability 
of  the  reconstruction  process  in  the  presence  of  noise  can  be  understood. 

For  this  model  of  signal  independent  noise  it  is  difficult  to  define  a 
signal  to  noise  ratio  for  a  given  noise  variance  and  Fresnel  zone  transform. 
As  mentioned  in  Chapt.  2,  as  the  Fresnel  zone  transform  moves  further 
away  from  the  Fourier  transform  plane  and  the  number  of  Fresnel  zones 
increases,  the  magnitude  data  tends  to  spread  out  while  the  phase 
information  takes  on  a  very  definite  circular  structure.  In  Sect.  5.2  it  was 
shown  that  the  probability  density  function,  and  thus  the  expected  value 
and  variance,  of  the  measured  signal  depends  on  the  ratio  of  the 
magnitude  of  the  true  signal  to  the  standard  deviation  of  the  noise.  For 
signal  independent  noise,  the  variance  of  the  noise  will  be  the  same  over 
the  entire  Fresnel  zone  transform,  but  the  magnitude  of  the  true  signal 
will  vary  from  pixel  to  pixel,  thus  the  signal  to  noise  ratio  will  vary  from 
pixel  to  pixel. 

The  values  of  the  maximum,  mean,  and  variance  of  the  magnitude 
information  for  several  different  Fresnel  zones  is  shown  in  table  5.1. 
Although  this  data  represents  the  Fresnel  zone  transforms  from  only  one 
image,  the  general  trends  of  a  decrease  in  the  maximum  value  and 
variance,  accompanied  by  an  increase  in  the  mean  value  has  been 
observed  for  other  images  that  were  examined.  It  should  also  be  noted 
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#  of  Fresnel 

Maximum 

Mean 

Variance 

Zones 

(xIOS) 

(x103) 

(x108) 

00.0 

64.5 

3.96 

8.54 

00.5 

37.0 

4.24 

8.52 

03.0 

5.57 

6.10 

8.33 

05.0 

3.78 

7.63 

8.12 

07.0 

2.80 

9.16 

7.86 

10.0 

2.13 

11.5 

7.38 

16.0 

1.45 

16.1 

6.12 

Table  5-1.  Values  for  the  maximum,  mean,  and  variance  of  the 
magnitude  of  the  Fresnel  zone  transform  for  several  different 
Fresnel  zones 

that  the  minimum  value  of  the  magnitude  data  is  generally  less  than  10, 
and  is  quite  often  equal  to  zero.  From  this  table  we  can  see  that  the 
known  magnitude  data  has  quite  a  large  dynamic  range,  and  thus  for  any 
particular  Fresnel  zone  transform,  the  signal  to  noise  ratio  will  vary 


significantly  over  the  transform  plane. 

To  investigate  the  effects  of  noise  on  reconstructing  from  Fresnel  zone 
magnitude-only  data,  various  amounts  of  noise  were  added  to  the  Fresnel 
zone  transform  at  several  different  Fresnel  zones.  Both  a  visual  inspection 


of  the  reconstructed  images  and  an  analytic  measure  of  the  error  in  the 
reconstruction  were  done. 


The  measure  of  the  absolute  error  in  the  transform  plane,  Em('),  is 
given  by 


1  N-l 
1 


0  n  =  0 


|G(m,n)|  -  |G.(m,n)| 


(5.3-1) 


a 
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where  i  is  the  iteration  number,  |G(m,n)j  is  the  measured  magnitude  data, 
and  |Gj(m,n)|  is  the  magnitude  data  calculated  from  the  Fresnel  zone 
transform  of  the  guess  for  the  image  after  i  iterations.  This  error  is  a  direct 
measure  of  how  well  the  iteration  is  progressing  in  terms  of  being 
consistent  with  the  known  magnitude  information. 

The  measure  of  the  relative  squared  error  in  the  image  plane,  EjO),  is 
given  by 


N-\  S-l 

I  I 

m—0  n=0 


g.im.n)  -  g.{m,n) 


S-l  N-1 

I  I 


2 


g^im.n) 


(5.3-2) 


fn  =  0  »i  =  0 

where  i  is  the  iteration  number,  gj(m,n)  is  the  calculated  guess  for  the 
image  after  i  iterations  and  before  the  image  plane  constraints  have  been 
imposed,  and  gi'(m,n)  is  the  result  of  applying  the  image  plane  constraints 
on  gi(m,n).  Although  this  error  is  very  difficult  to  relate  to  the  known 
transform  magnitude  information,  this  error  does  give  an  indication  as  to 
how  much  the  guess  for  the  image  is  having  to  be  modified  to  fulfill  the 
image  plane  constraints. 

For  the  case  when  there  is  no  noise  in  the  measured  data,  it  is  possible 
to  fulfill  all  the  constraints  imposed  in  both  the  image  and  transform 
domains,  so  that  these  errors  can  go  to  zero.  When  noise  is  added  to  the 
measured  data,  the  magnitude  information  imposed  in  the  transform 
domain  may  no  longer  be  compatible  with  the  constraints  of  a  real 
valued,  non-negative  object.  The  absolute  error  in  the  transform  domain 
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and  the  relative  squared  error  in  the  image  plane  could  not  both  go  to 
zero  in  this  case. 

Figure  5-8  is  a  plot  of  the  absolute  error,  Em('),  for  the  first  75  iterations 
at  7.0  Fresnel  zones  for  various  amounts  of  noise.  We  see  in  these  plots 
that  the  error  tends  to  change  the  most  in  the  first  several  iterations  and 
then  decreases  very  slowly.  Note  also  that  the  rate  of  change  of  the  error 
from  iteration  2  until  approximately  iteration  25  is  directly  proportional  to 
the  amount  of  noise  in  the  data.  For  approximately  the  first  25  iterations, 
the  more  noise  there  is  in  the  measured  data,  the  greater  decrease  there  is 
in  the  absolute  error  between  the  measured  noisy  data  and  the  calculated 
transform  magnitude.  Although  Fig.  5-8  contains  the  error  plots  for  only 
one  realization  of  the  noise  at  7.0  Fresnel  zones,  the  general  structure  and 
shape  of  the  error  curves  are  very  similar  both  for  different  realizations  of 
the  noise  and  for  data  recorded  at  different  locations  in  the  Fresnel  zone. 

One  of  the  observations  that  is  made  in  examining  the  plots  shown  in 
Fig.  5-8  is  that  the  error  tends  to  stagnate  on  a  level  which  is 
approximately  equal  to  one-half  the  value  of  the  standard  deviation  of 
the  signal  independent  noise.  Figure  5-9  is  a  plot  of  the  value  of  this 
stagnation  error  verses  the  standard  deviation  of  the  noise  for  several 
different  Fresnel  zones.  The  value  of  the  error  at  the  stagnation  level  is 
calculated  as  the  average  value  of  the  error  between  iteration  50  and 
iteration  75.  From  this  plot  we  can  see  that  the  error  level  is 
approximately  equal  to  one-half  the  value  of  the  standard  deviation  of 
the  noise,  independent  of  the  Fresnel  zone  at  which  the  data  was 
measured. 
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Fig.  5-8  Absolute  error  between  current  transform  magnitude  and  measured  transform 
magnitude  for  magnitude-only  reconstruction  at  7.0  Fresnel  zones.  Errors  are 
shown  for  cases  when  the  standard  deviation  of  the  noise  is  equal  to  0,  5000, 
10000, 15000, 20000,  25000,  and  30000. 

Figure  5-10  is  a  plot  of  the  relative  squared  error,  Ei(i),  for  the  first  75 
iterations  at  7.0  Fresnel  zones  for  various  amounts  of  noise.  We  see  in 
these  plots  that  the  error  tends  to  change  the  most  in  the  first  several 
iterations  and  then  decreases  very  slowly.  Note  also  that  the  rate  of 
change  of  the  error  for  the  first  few  iterations  is  directly  proportional  to 
the  amount  of  noise  in  the  data.  For  approximately  the  first  25  iterations, 
the  more  noise  there  is  in  the  measured  data,  the  greater  decrease  there  is 
in  the  relative  squared  error  in  the  image  plane. 
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Fig.  5-9  Average  value  of  the  absolute  error  in  the  magnitude-only  reconstruction 
for  iterations  50  through  75  for  various  amounts  of  noise  at  several 
different  Fresnel  zone  locations.  The  plotted  symbols  1-6  are  for  Fresnel 
zones  equal  to  0.0, 0.5, 3.0, 5.0, 7.0,  and  10.0  respectively. 

As  shown  in  Fig.  5-10,  the  error  tends  to  stagnate  at  a  level  which  is 
directly  proportional  to  the  value  of  the  standard  deviation  of  the  signal 
independent  noise.  Figure  5-1 1  is  a  plot  of  the  value  of  this  steady  error 
verses  the  standard  deviation  of  the  noise  for  several  different  Fresnel 
zones.  The  value  of  the  error  at  the  stagnation  level  is  calculated  as  the 
average  value  of  the  error  between  iteration  50  and  iteration  75.  From 
this  plot  we  can  see  that  the  stagnation  error  level  for  a  given  noise 
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Fig.  5-10  Relative  squared  error  in  the  image  plane  for  75  iterations  at  7.0  Fresnel 
zones  for  magnitude-only  reconstructions.  Errors  are  shown  for  the  cases 
when  the  standard  deviation  of  the  noise,  o,  is  equal  to  0,  5000, 1 0000, 1 5000, 
20000, 25000,  and  30000. 

variance  is  independent  of  the  Fresnel  zone  at  which  the  data  was 
measured. 

To  see  how  this  absolute  error  in  the  reconstruction  relates  to  the 
quality  of  the  reconstructed  image,  consider  the  sample  of  images  shown 
in  Figs.  5-12  through  5-14.  Figure  5-12  is  an  example  of  a  reconstruction  at 
3.0  Fresnel  zones  after  15  and  75  iterations  from  data  corrupted  by  noise 
with  standard  deviations  of  0,  2500,  10000,  and  25000.  From  these 
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Fig.  5-11  Average  value  of  the  relative  squared  error  in  the  magnitude-only 
reconstruction  for  iterations  50  through  75  for  various  amounts  of  noise  at 
several  different  Fresnel  zone  locations.  The  plotted  symbols  1-6  are  for 
Fresnel  zones  0.0, 0.5, 3.0, 5.0, 7.0,  and  10.0  respectively. 

reconstructed  images  we  see  that  the  quality  of  the  reconstruction 
changes  less  and  less  from  iteration  15  to  iteration  75  as  the  amount  of 
noise  is  increased.  For  noise  values  equal  to  0  and  2500  there  is  a  visible 
increase  in  the  quality  of  the  reconstruction,  for  noise  with  a  standard 
deviation  of  10000  there  is  no  visible  change  in  this  barely  recognizable 
image,  and  for  noise  with  a  standard  deviation  of  25000  the  image  is  not 
recognizable  for  any  number  of  iterations. 
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Fig.  5-12  Magnitude-only  reconstructions  from  noisy  data  at  3  0  Fresnel  zones  The 


reconstructions  are  shown  after  15  iterations  (left)  and  75  iterations  (right) 
for  various  amounts  of  noise  as  indicated 


OS  10000 


o  s  25000 


Fig.  5-13  Magnitude-only  reo 
reconstructions  are 
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Fig.  5-14  Magnitude-only  reconstructions  from  noisy  data  at  10.0  Fresnel  zones.  Th* 
reconstructions  are  shown  after  15  iterations  (left)  and  75  iterations  (right) 
for  various  amounts  of  noise  as  indicated. 


5.3.2  Phase-only  Reconstruction 

The  noise  model  developed  in  Sect.  5.2  was  also  applied  to  various 
Fresnel  zone  transforms  before  calculating  the  phase  angle  data,  and  then 
attempting  to  reconstruct  an  image  from  this  noisy  data. 

As  discussed  earlier  in  Sect.  5.3.1,  it  is  difficult  to  define  a  signal  to 
noise  ratio  for  a  given  noise  variance  and  Fresnel  zone  transform,  in  Sect. 

5.2  it  was  shown  that  the  probability  density  function,  and  thus  the 
expected  value  and  variance  of  the  measured  phase  angle  data  depends 
on  the  ratio  of  the  magnitude  of  the  true  signal  to  the  standard  deviation 
of  the  noise,  in  phase-only  reconstruction  it  is  assumed  that  this 
magnitude  information  is  not  known.  In  addition,  the  calculation  of  the 
phase  information  in  the  Fresnel  zone  transform  may  require  the 
measurement  of  several  intensity  distributions  or  interference  patterns, 
followed  by  the  numerical  calculation  of  the  phase  angle  at  each  pixel. 
Although  the  measurement  of  the  phase-only  information  is  often  a 
complicated  one,  the  noise  model  chosen  has  the  property  that  the 
stronger  the  signal  at  a  given  pixel,  the  less  corrupted  the  phase 
information  is  at  that  pixel  for  a  given  amount  of  noise. 

To  investigate  the  effects  of  noise  on  reconstructing  an  image  from 
Fresnel  zone  phase-only  data,  various  amounts  of  noise  were  added  to  the 
Fresnel  zone  transform  at  several  different  Fresnel  zones.  Both  a  visual 
inspection  of  the  reconstructed  images  and  an  analytic  measure  of  the 
error  in  the  reconstruction  were  done. 

A  measure  of  the  analytic  error  in  the  transform  domain,  Ep(i),  is 
defined  as  by 


exp{i  0(m,n)}  —  exp{i  0 .  (m,n)} 


(5.3-3) 


S-l  N-l 
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where  exp{iO(m,n)}  is  the  knowr)  phase  information,  and  exp{i0i(m,n)}  is 
the  phase  of  the  Fresnel  zone  transform  of  the  guess  for  the  image  after  i 
iterations.  The  relative  squared  error  in  the  image  plane,  defined  in  Eq. 
(5.3-2),  was  used  as  an  analytic  measure  of  the  error  in  the  image  plane 
during  the  reconstruction  process. 

A  plot  of  the  absolute  phase  error,  Ep(i).  is  shown  in  Fig.  5-15.  Although 
the  presence  of  only  a  small  amount  of  noise  results  in  a  large  increase  in 
the  absolute  phase  error,  there  is  not  a  large  degradation  in  the  visual 
quality  of  the  reconstructed  image.  It  is  interesting  to  note  that  when 
noise  is  present  in  the  known  data,  the  absolute  phase  error  has  a 
minimum  after  3  or  4  iterations,  then  increases  to  a  steady  state  value 
which  is  slightly  greater  than  the  minimum  value  located  near  iteration 
number  3. 

A  plot  of  the  stagnation  value  of  the  absolute  phase  error  verses  the 
standard  deviation  of  the  noise  is  shown  in  Fig.  5-16.  The  stagnation  value 
of  the  absolute  phase  error  is  equal  to  the  average  value  of  the  error  for 
iterations  15  through  20.  From  Fig.  5-16  we  see  that  even  a  small  amount 
of  noise  results  in  a  large  increase  in  the  absolute  phase  error,  and  that  the 
phase  error  is  slightly  less  at  locations  which  are  farther  away  from  the 
Fourier  transform  plane. 

Figure  5-17  is  a  plot  of  the  relative  squared  error,  Ei(i),  for  the  first  20 
iterations  at  7.0  Fresnel  zones  for  various  amounts  of  noise.  We  see  from 
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Fig.  5-15  Absolute  phase  error  in  the  transform  plane  for  20  iterations  at  7.0  Fresnel 
zones  from  phase-only  data.  Errors  are  shown  for  cases  where  the  standard 
deviation  of  the  noise,  o,  is  equal  to  0,  2500, 15000,  and  30000. 

these  plots  that  the  error  tends  to  change  the  most  in  the  first  few 
iterations  and  then  decrease  very  slowly.  The  relative  squared  error  is  not 
very  sensitive  to  the  preserice  of  small  amounts  of  noise. 

Figure  5-18  is  a  plot  of  the  average  value  of  the  relative  squared  error 
from  iteration  number  15  to  iteration  number  20  verses  the  standard 
deviation  of  the  noise  for  several  different  Fresnel  zones.  From  this  plot 
we  see  that  the  stagnation  value  of  the  relative  squared  error  increases 
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Fig.  5-16  Average  value  of  the  absolute  phase  error  in  the  image  plane  for  iterations 
15  through  20  at  Fresnel  zones  equal  to  0.0,  0.5,  3.0,  5.0,  7.0,  and  10.0 
(indicated  by  the  symbols  1-6  respectively).  Errors  are  plotted  for  cases  where 
the  standard  deviation  of  the  noise,  o,  is  equal  to  0,  2500,  5000,  7500,  10000, 
15000, 20000, 25000,  and  30000. 

both  as  the  amount  of  noise  increases,  and  as  the  number  of  Fresnel  zones 
increases. 

To  see  how  this  error  in  the  reconstruction  relates  to  the  quality  of  the 
reconstructed  image,  consider  the  sample  of  images  shown  in  Figs.  5-19 
through  5-21 .  Figure  5-19  is  an  example  of  a  reconstruction  at  3.0  Fresnel 
zones  after  3  and  20  iterations  from  data  corrupted  by  noise  with 
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Fig.  5-17  Relative  squared  error  in  the  image  plane  for  75  iterations  at  7.0  Fresnel 
zones  from  phase-only  data.  Errors  are  shown  for  cases  where  the  standard 
deviation  of  the  noise,  o,  is  equal  to  0, 2500, 15000,  and  30000. 

deviations  of  0, 2500, 1 5000,  and  30000.  From  these  reconstructed  images 
we  see  that  the  quality  of  the  reconstruction  changes  both  as  the  number 
of  iterations  increases,  and  as  the  amount  of  noise  is  increased.  Generally, 
in  the  first  few  iterations  the  fine  detail  and  sharp  edges  show  up  with  a 
higher  contrast  than  there  is  in  iteration  numbers  10  through  20.  The 
reconstruction  procedure  is  forfeiting  some  of  the  high  contrast  detail  in 
an  attempt  to  remove  the  low  frequency  noise  ripples.  This  smoothing  of 
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Fig.  5-18  Stagnation  level  of  the  relative  squared  error  in  phase-only  reconstruction  for 
several  amounts  of  noise  at  different  Fresnel  zone  locations.  The  plotted 
symbols  1-6  are  for  Fresnel  zones  equal  to  0.0,  0.5,  3.0,  5.0,  7.0,  andlO.O 
respectively. 


the  high  contrast  detail  to  satisfy  the  image  plane  constraints  results  in  a 
slight  increase  in  the  absolute  phase  error  that  is  shown  in  Fig.  5-15.  The 
end  result  is  a  blurry  but  easily  recognized  image.  Even  in  the  presence  of 
large  amount  of  noise,  reconstruction  from  phase-only  information  yields 
a  surprisingly  good  reconstructed  image. 

Figures  5-20  and  5-21  show  the  same  set  of  images  discussed  for  Fig.  5- 
19,  but  the  reconstructions  for  Fig.  5-20  were  done  at  7.0  Fresnel  zones. 
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and  the  reconstructions  for  Fig.  5-21  were  done  at  10.0  Fresnel  zones. 
From  these  images  we  observe  that  as  the  number  of  Fresnel  zones 
increases  the  reconstruction  algorithm  is  very  stable,  and  indeed  is 
insensitive  to  the  presence  of  even  large  amount  of  noise. 
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5.4  Coherent  Multiplicative.  Signal-dependent  Noise 
5.4.1.  Noise  Model 

In  Sect.  5.2  a  model  for  coherently  additive,  signal  independent, 
uncorrelated  Gaussian  noise  was  developed.  In  this  model,  the  variance  of 
the  noise  was  independent  of  the  magnitude  of  the  signal  at  each  pixel, 
and  the  variance  of  the  noise  was  a  constant  over  the  entire  Fresnel  zone 
transform.  For  the  case  of  signal  dependent,  or  multiplicative,  noise,  the 
variance  of  the  noise  will  depend  upon  the  magnitude  of  the  signal,  and 
thus  the  variance  of  the  noise  will  vary  across  the  Fresnel  zone  transform. 

In  this  multiplicative  noise  model,  a  noise  signal  is  multiplied  by  the 
magnitude  of  the  true  signal,  and  then  added  coherently  to  the  true 
signal.  The  additive  noise  signal  will  be  represented  as 

n  =r|u>J(x  +  or) ,  (5.4-1) 

where  x  and  y  are  independent,  uncorrelated,  Gaussian  distributed, 
random  variable  with  zero  mean  and  unit  variance,  |ws|  is  the  magnitude 
of  the  true  signal,  and  r  is  a  scaling  factor  used  to  control  the  signal  to 
noise  ratio. 

Recording  devices  such  as  film  or  photo-detectors  respond  to  the 
intensity  of  the  signal,  |ws|2,  and  the  film  grain  noise  or  shot  noise  present 
in  the  detection  system  will  have  a  standard  deviation  which  is 
proportional  to  the  square  root  of  the  detected  intensity.  Since  x  and  y 
come  from  a  distribution  with  unit  variance,  multiplying  by  |ws|  will  give  a 
noise  value  with  a  standard  deviation  equal  to  the  square^  root  of  the 
intensity.  Multiplying  the  noise  by  the  variable  r  will  provide  a  means  of 
choosing  the  signal  to  noise  ratio.  If  the  signal  intensity  is  |wsp  and  the 
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variance  of  the  noise  is  r‘|ws|2  .then  the  signal  to  noise  ratio  will  be 
defined  as 

i  =  J-  (5.4-2) 

N  r* 

The  graphical  model  and  analytic  derivations  presented  in  Sect.  5.2  for 
the  case  of  signal  independent  additive  noise  are  still  valid  for  the  case  of 
signal  dependent  multiplicative  noise.  The  only  difference  is  that  for 
signal  independent  noise,  the  standard  deviation  of  the  noise,  o,  is  an 
independent  variable,  and  now  for  the  case  of  signal  dependent  noise, 
the  standard  deviation  of  the  noise  is  given  by: 

o  =  r|u>J,  (5.4-3) 

where  r  is  used  to  vary  the  amount  of  noise.  The  standard  deviation  of 
the  noise  varies  from  pixel  to  pixel  in  the  Fresnel  zone  transform  and  is 
equal  to  the  product  of  r  and  the  magnitude  of  the  signal  at  the  given 
pixel. 

The  joint  probability  density  function  for  the  magnitude  and  phase  will 
be 


2  n(r  s) ' 


exp  -[(6co8e -8)‘  +  (6  8in0)'l  /ZCFb)' 


(5.4-4) 


6  2  0  and,  —a  <  0  S  re  . 


The  probability  density  function  for  the  measured  magnitude,  b,  is  given 
by 


A(6)  = 


(TsV 


exp  j  -(6^  +  s^)  /2(rs)^ 


62  0  , 


(5.4-5) 
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and  the  probability  density  function  for  the  measured  phase  angle,  6,  is 
given  by 

|,  +  2 .^2^)1  (5.4-6) 

— n  <  0  S  n  , 

The  plots  for  the  probability  density  functions  of  the  magnitude  and 
phase  are  the  same  as  the  plots  shown  in  Figs.  5-4  and  5-6  with  the 
understanding  that  o^sr  and  the  parameter  ks  s/a  is  to  be  replaced  by 
k  =  l/r. 


5.5  Effects  of  Multiplicative  Noise  on  Image  Reconstruction 
5.5.1 .  Magnitude-only  Reconstruction 

For  the  model  discussed  in  Sect.  5.4.  a  signal  to  noise  ratio  (S/N) 
denoted  by 

S/N  =  lfT^  ,  (5.5-1) 

where  F  is  defined  in  Eq.  (5.4-1)  as  a  method  of  controlling  the  signal  to 
noise  ratio,  has  been  used.  To  investigate  the  effects  of  noise  on 
reconstructing  an  image  from  Fresnel  zone  magnitude-only  data,  various 
amounts  of  noise  were  added  to  the  Fresnel  zone  transform  at  several 
different  locations.  Both  a  visual  inspection  of  the  reconstructed  images 
and  an  analytic  measure  of  the  error  in  the  reconstruction  were  done. 

A  measure  of  how  well  the  reconstruction  is  progressing  in  terms  of 
being  consistent  with  the  known  magnitude  information  is  the  absolute 
error  in  the  transform  plane,  Eivi('^  defined  in  Eq.(5.3-1).  Figure  5-22  is  a 
plot  of  the  absolute  error,  Em(')»  for  the  first  75  iterations  at  7.0  Fresnel 
zones  for  various  amounts  of  noise.  We  can  see  from  these  plots  that  the 
error  drops  sharply  in  the  first  two  iterations,  and  then  stagnates  at  a 
value  which  is  proportional  to  the  amount  of  noise  present  in  the  data. 
For  the  case  when  no  noise  is  present,  the  error  continues  to  decrease 
steadily  toward  zero. 

The  value  at  which  the  error  levels  off  depends  upon  both  the  number 
of  Fresnel  zones  and  the  signal  to  noise  ratio.  Figure  5-23  is  a  plot  of  the 
average  value  of  the  absolute  error  for  iterations  50  through  75  for 
various  amounts  of  noise  at  several  different  Fresnel  zone  locations.  The 
plotted  symbols  1-6  correspond  to  the  number  of  Fresnel  zones  being 
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Fig.  S-22  Absolute  error  between  current  transform  magnitude  and  measured 
transform  magnitude  for  magnitude-only  reconstruction  at  7.0  Fresnel  zones. 
Errors  are  shown  for  cases  when  the  signal  to  noise  ratio  is  equal  to  50.0, 
10.0, 5.0, 1.0, 0.5,  and  0.1. 

equal  to  0.0,  0.5,  3.0,  5.0,  7.0,  and  10.0  respectively.  These  errors  are 
plotted  for  various  values  of  the  parameter  F  corresponding  to  signal  to 
noise  ratios  of  «,  50.0, 10.0, 5.0, 1 .0, 0.5,and  0.1.  From  this  plot  we  see  that 
the  error  level  depends  both  on  the  signal  to  noise  ratio  and  on  the 
number  of  Fresnel  zones.  The  error  level  increases  both  with  the  amount 
of  noise  (i.e.  a  decrease  in  S/N)  and  with  the  number  of  Fresnel  zones. 
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The  error  levels  plotted  in  Fig.  5-23  are  close  to  the  value  of  r/2  times 
the  average  value  of  the  true  transform  magnitude  data.  The  average 
standard  deviation  of  the  noise  in  this  signal  dependent  noise  model  is  T 
times  the  average  value  of  the  transform  magnitude  data.  This  error  level 
corresponds  very  well  with  the  case  of  signal  independent  noise,  where 
the  error  level  was  found  to  be  approximately  equal  to  o/2. 


Fig.  5-23  Average  values  of  the  absolute  error  in  the  magnitude-only  reconstruction 
for  iterations  50  through  75  at  Fresnel  zones  equal  to  0.0,  0.5,  3.0,  5.0,  7.0, 
and  10.0  (shown  by  the  characters  1-6  respectivley).  Errors  are  shown  for  the 
cases  when  the  signal  to  noise  ratio  is  equal  to  »,  50.0,  10.0,  3.0,  1 .0,  0.5,  and 
0.1. 
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For  completeness  and  for  comparison  with  the  case  of  signal 
independent  noise,  a  plot  of  the  relative  squared  error  in  the  image  plane 
is  shown  in  Fig.  5-24.  The  basic  structure  in  these  plots;  where  the  error 
decreases  most  in  the  first  several  iterations,  and  then  levels  off  to  a  value 
which  is  proportional  to  the  signal  to  noise  ratio,  is  common  to  both  other 
realizations  of  the  noise,  and  to  reconstructions  at  different  Fresnel  zone 
locations. 

A  plot  of  the  stagnation  value  of  the  relative  squared  error  in  the 
image  plane  is  shown  in  Fig.  5-25.  For  a  noise  parameter,  P,  which  is  less 
than  or  equal  to  one  (i.e.  S/N  ^1.0)  the  relative  squared  error  is 
independent  of  the  number  of  Fresnel  zones.  For  large  amounts  of  noise 
(i.e.  r>  1 .),  the  error  depends  both  on  the  number  of  Fresnel  zones  and  on 
the  signal  to  noise  ratio.  The  plots  shown  in  Fig,  5-25  represent  only  one 
realization  for  the  noise.  Different  choices  for  the  noise  yield  error  plots 
which  are  very  similar  to  Fig,  5-25  for  small  noise  levels,  but  for  large 
amount  of  noise  the  stagnation  value  of  the  noise  can  change  by 
approximately  ±0.1. 

To  see  how  the  analytic  error  relates  to  the  visual  quality  of  the 
reconstructed  image,  consider  the  sample  of  images  shown  in  Figs.  5-26 
through  5-28.  Figure  5-26  is  an  example  of  a  reconstruction  at  3.0  Fresnel 
zones  after  1 5  and  75  iterations  from  data  with  a  signal  to  noise  ratio  of  «, 
50.0,  5.0,  and  0.1.  From  these  reconstructed  images  we  see  that  the 
recorstruction  does  not  change  noticeably  from  iteration  15  to  iteration 
75,  which  is  consistent  with  the  error  plots  shown  in  Figs.  5-22  and  5-24. 
With  a  signal  to  noise  ratio  of  50.0,  the  reconstructed  image  is  only  slightly 
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Fig.  5-24  Relative  squared  error  in  the  image  plane  for  75  iterations  at  7.0  Fresnel 
zones  for  magnitude-only  reconstuctions.  Errors  are  shown  for  the  cases 
when  the  signal  to  noise  ratio  is  equal  to®,  50.0, 10.0,  3.0, 1.0, 0.5,  and  0.1. 

degraded.  With  a  signal  to  noise  ratio  of  5.0,  the  reconstructed  image  is 
barely  recognizable.  With  a  signal  to  noise  ratio  of  0.1,  the  image  is 
essentially  unrecognizable. 

Figures  5-27  and  5-28  are  the  same  set  of  reconstructed  images  as 
discussed  for  Fig.  5-26,  but  the  reconstructions  for  Fig.  5-27  were  done  at 
7.0  Fresnel  zones,  and  the  reconstructions  in  Fig.  5-28  were  done  at  10.0 
Fresnel  zones. 


i| 
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Fig.  5-25  Average  values  of  the  relative  squared  error  in  the  image  plane  for  iterations 
50  through  75  at  Fresnel  zones  equal  to  0.0,  0.5,  3.0,  5.0,  7.0,  and  10.0  (shown 
by  the  characters  1-6  respectivley).  Errors  are  shown  for  the  cases  when  the 
signal  to  noise  ratio  is  equal  to  50.0, 1 0.0, 3.0, 1 .0,  0.5,  and  0, 1 


5.5.2.  Phase-only  Reconstruction 

To  investigate  the  effects  of  noise  on  reconstructing  an  image  from 
Fresnel  zone  phase-only  data,  various  amounts  of  noise  were  added  to  the 
Fresnel  zone  transform  at  several  different  locations.  Both  a  visual 
inspection  of  the  reconstructed  images  and  an  analytic  measure  of  the 
error  in  the  reconstruction  were  done. 

The  absolute  phase  error  in  the  transform  plane,  Ep(>),  defined  in  Eq. 
(5.3-3)  and  the  relative  squared  error  in  the  image  plane,  EjO),  defined  in 
Eq.  (5.3-2)  were  used  as  analytic  measures  of  the  errors  during  the 
reconstruction  procedure. 

Figure  5-29  is  a  plot  of  the  absolute  phase  error,  EpO),  for  the  first  20 
iterations  at  7.0  Fresnel  zones  for  various  signal  to  noise  ratios.  The  error 
at  the  first  iteration  is  approximately  .5,  and  is  independent  of  the  amount 
of  noise.  The  decrease  in  the  error  during  the  first  few  iterations  is 
indirectly  proportional  to  the  amount  of  noise  that  is  present,  when  no 
noise  is  present  the  error  continues  to  decrease  even  after  20  iterations. 
For  the  cases  with  large  amounts  of  noise  the  error  decreases  slightly  after 
the  first  few  iterations  and  then  stagnates  at  a  constant  value. 

A  plot  of  the  stagnation  values  of  the  absolute  phase  error  for  several 
different  Fresnel  zones  is  shown  in  Fig.  5-30.  From  this  plot  we  see  that  for 
small  amounts  of  noise,  the  error  level  is  directly  proportional  to  the 
number  of  Fresnel  zones.  When  the  signal  to  noise  ratio  is  less  than  or 
equal  to  one  (i.e.  1.)  the  error  level  is  independent  of  the  number  of 

Fresnel  zones,  and  in  all  cases  is  very  close  to  0.5.  The  error  in 
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Fig.  5-29  Absolute  phase  error  in  the  transform  plane  for  20  iterations  at  7.0  Fresnel 
zones  from  phase-only  data.  Errors  are  shown  for  cases  where  the  standard 
deviation  of  the  noise,  o,  is  equal  to  0. 2500, 1 5000,  and  30000. 

reconstructing  from  Fourier  transform  phase-only  data  is  consistently 
lower  than  the  error  in  reconstructing  from  Fresnel  zone  phase-only  data. 

Figure  5-31  is  a  plot  of  the  relative  squared  error,  Ei(i),  for  the  first  20 
iterations  at  7.0  Fresnel  zones  for  various  signal  to  noise  ratios,  in  the 
presence  of  noise  the  error  decreases  for  the  first  few  iterations,  and  then 
stagnates  at  a  level  which  is  indirectly  proportional  to  the  signal  to  noise 
ratio. 
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Fig.  5-30  Average  values  of  the  absolute  phase  error  for  iterations  15  through  20  at 
Fresnel  zones  equal  to  0.0. 0.5.  3.0.  5.0. 7.0.  and  10.0  (shown  by  the  characters 
1-6  respectivley).  Errors  are  shown  for  the  cases  when  the  signal  to  noise  ratio 
is  equal  to  50.0. 1 0.0. 3.0. 1 .0. 0.5.  and  0. 1 . 


A  plot  of  the  stagnation  values  of  the  relative  squared  error  for  several 
different  Fresnel  zones  is  shown  in  Fig.  5-32.  For  the  cases  when  the  signal 
to  noise  parameter,  F,  is  greater  than  or  equal  to  .5,  the  error  is  slightly 
less  in  the  Fourier  transform  plane.  When  the  parameter  F  is  less  than  .5, 
the  stagnation  level  of  the  error  is  independent  of  the  number  of  Fresnel 
zones. 
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Fig.  5-31  Relative  squared  error  in  the  image  plane  for  20  iterations  at  7.0  Fresnel 
zones  for  phase-only  reconstructions.  Errors  are  shown  for  the  cases  when 
the  signal  to  noise  ratio  is  equal  to  50.0,  5.0,  and  0. 1 . 

To  see  how  this  error  in  the  reconstruction  relates  to  the  quality  of  the 
reconstructed  image,  consider  the  sample  of  images  shown  in  Figs.  5-33 
through  5-35.  Figure  5-33  is  an  example  of  a  reconstruction  at  3.0  Fresnel 
zones  after  3  and  20  iterations  from  data  corrupted  by  noise  with  signal 
to  noise  ratio  of  »,  50.0,  5.0,  and  0.1..  From  these  reconstructed  images  we 
see  that  the  quality  of  the  reconstruction  changes  both  as  the  number  of 
iterations  increases,  and  as  the  amount  of  noise  is  increased.  Generally,  in 
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Fig.  5-32  Average  value  of  the  relative  squared  error  in  the  image  plane  for  iterations 
15  through  20  at  Fresnel  zones  equal  to  0.0,  0.5,  3.0,  5.0,  7.0,  and  10.0 
(indicated  by  the  characters  1-6  respectively).  Errors  are  plotted  for  the  cases 
when  the  sigani  to  noise  ratio  is  equal  to  «,  50.0,  5.0,  and  0. 1 . 


the  first  few  iterations  the  fine  detail  and  sharp  edges  show  up  with  a 
higher  contrast  than  there  is  in  iteration  numbers  10  through  20.  The 
reconstruction  procedure  is  forfeiting  some  of  the  high  contrast  detail  in 
an  attempt  to  remove  the  low  frequency  noise  ripples. 

Figures  5-34  and  5-35  are  the  same  set  of  images  discussed  for  Fig.  5-33, 
but  the  reconstructions  for  Fig.  5-34  were  done  at  7.0  Fresnel  zones,  and 
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the  reconstructions  for  Fig.  5-35  were  done  at  10.0  Fresnel  zones.  From 
these  images  we  observe  that  as  the  number  of  Fresnel  zones  increases  the 
reconstruction  algorithm  is  very  stable.and  indeed  is  very  insensitive  to  the 
presence  of  even  large  amount  of  noise. 
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Fig.  5-34  Phase-only  reconstructions  from  noisy  data  at  7.0  Fresnel  zones.  The 
reconstructions  are  shown  after  3  iterations  (left)  and  20  iterations  (right)  for 
various  signal  to  noise  ratios  as  indicated. 
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Fig.  5-35  Phase-only  reconstructions  from  noisy  data  at  10.0  Fresnel  zones.  The 
reconstructions  are  shown  after  3  iterations  (left)  and  20  iterations  (right)  for 
various  signal  to  noise  ratios  as  indicated. 
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Chapter  Six 


Effects  of  an  Error  in  the  Number  of  Fresnel  Zones 

6.1  Introduction 

Up  to  this  point,  the  reconstruction  algorithms  have  assumed  that  the 
number  of  Fresnel  zones,  Zf,  is  a  known  quantity.  The  number  of  Fresnel 
zones  is  connected  to  the  distance  out  of  the  Fourier  transform  plane  by 
the  relation 


where  di  is  the  distance  behind  the  lens  (see  Fig.  2-1),  F  is  the  focal  length 
of  the  lens,  A  is  the  wavelength,  N  is  the  number  of  sample  points,  and  T  is 
the  sampling  interval.  From  this  equation  we  can  see  that  if  the  distance 
di  is  given  incorrectly,  the  result  is  an  error  in  the  parameter  Zp.  An  error 
similar  to  this  was  also  discussed  by  Childs  and  Misell  [1973]  for  the  case  of 
an  error  in  the  amount  of  defocus  in  image  deconvolution,  where  it  was 
discovered  that  "the  use  of  the  correct  defocus  value  is  extremely 
important." 

An  incorrect  value  of  di  could  be  the  result  of  experimental  error, 
where  the  distance  is  either  measured  incorrectly,  or  can  not  be  measured 
accurately.  In  these  cases,  one  should  be  able  to  estimate  a  value  for  Zf, 
but  would  like  to  know  the  effects  of  having  to  use  an  estimated  value. 
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Section  6.2  will  examine  the  effects  upon  the  reconstructed  images  of  this 
type  of  error  in  the  reconstruction  process. 

A  different  problem  exists  when  the  given  information,  either 
magnitude-only  or  phase-only  data,  is  from  an  unknown  location  in  the 
Fresnel  zone  region.  This  may  be  the  case  when  the  distance  di  is  simply 
unavailable  for  measurement.  An  alternate  possibility  is  that  the  data  has 
been  purposely  encrypted  by  a  method  similar  to  the  technique  suggested 
by  VanderLught  [1985].  Section  6.3  will  describe  a  method  of  searching 
through  the  Fresnel  zone  region,  and  thus  reconstruct  an  image  from 
partial  data  at  an  unknown  location  in  the  Fresnel  region. 
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6.2  Use  of  Incorrect  Number  of  Fresnel  Zones 

As  discussed  in  Chapter  2,  the  reconstruction  algorithms  require  the 
use  of  a  complex  Gaussian  phase  term,  exp{-inai(x2  +  y2)}.  This  complex 
Gaussian  is  not  needed  for  an  inverse  Fresnel  zone  transform  if  the  object 
is  known  to  be  real  and  non-negative,  but  this  term  is  crucial  to  the 
calculation  of  the  forward  Fresnel  zone  transform  of  an  image.  This 
section  will  investigate  the  effects  of  using  an  estimate  for  the  number  of 
Fresnel  zones,  that  is,  the  value  of  If  used  in  the  reconstruction  algorithm 
is  nearly  equal  to  the  value  of  If  which  gave  rise  to  the  given  partial 
information. 

Figure  6-1  is  the  original  256  X  256  X  8-bit  image.  The  corresponding 


Fig.  6-1 .  The  original  digitized  image,  256  X  256  X  8-bits. 

magnitudes  of  the  Fresnel  zone  transforms  at  7.3,  7.7,  and  8.1  Fresnel 
zones  are  shown  in  Figs.  6-2(a),  (b),  and  (c)  respectively.  The  phases  of  the 
Fresnel  zone  transforms  at  7.3,  7.7,  and  8.1  Fresnel  zones  are  shown  in 
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(a)ZF  «  7.3 


I 

(b)ZF  =  7.7 
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{c)Zf  -  8.1 
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Fig  6-2. 
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Magnitude  of  the  Fresnel  zone  transforms  at  (a)7.3,  (b)7.7,  and  (c)8. 1  Fresnel 
zones 
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Figs.  6-3(a),  (b),  and  (c)  respectively.  The  Fresnel  zone  values  of  7.3  and  8.1 
correspond  to  differences  of  approximately  5%  from  the  central  value  of 
7.7  Fresnel  zones.  From  Figs.  2-4  and  2-9  it  is  seen  that  for  a  given  image, 
both  the  phase  and  magnitude  of  the  Fresnel  zone  transform  depend 
upon  the  number  of  Fresnel  zones.  The  change  in  the  phase  and 
magnitude  as  one  moves  through  the  Fresnel  zone  region  is  a  smooth  and 
gradual  transition.  While  there  is  a  visually  discernable  difference  in  the 
Fresnel  zone  transforms  from  locations  far  apart  in  the  Fresnel  zone 
region,  the  change  is  slight  for  locations  near  one  another  in  the  Fresnel 
zone  region.  Note  that  the  magnitudes  and  phases  shown  in  Figs.  6-2  and 
6-3  differ  only  slightly  for  the  different  values  of  the  number  of  Fresnel 
zones. 

6.2.2  Phase-only  Reconstruction 

To  examine  the  sensitivity  of  the  reconstruction  algorithm,  three 
reconstructions  using  phase-only  data  from  7.7  Fresnel  zones  were 
performed.  The  parameter  Zf  was  set  equal  to  7.3,  7.7  and  8. 1  for  each  of 
the  attempted  reconstructions.  The  values  of  7.3  and  8.1  correspond  to 
the  case  when  data  was  actually  obtained  at  7.7  Fresnel  zones,  but  the 
data  is  incorrectly  assumed  to  be  the  result  of  a  Fresnel  zone  transform 
which  differs  by  5%  from  the  correct  value.  Figures  6-4(a),  (b),  and  (c)  are 
the  results  of  these  phase-only  reconstructions  after  10  iterations.  Note 
that  there  is  still  a  recognizable  image  reconstructed  but  the  quality  of  the 
reconstruction  is  below  that  which  is  obtained  when  the  correct  number 
of  Fresnel  zones  is  used. 
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(a)ZF  «  7.3 


(b)ZF  *  7.7 


(c)Zf  =  8.1 


Fig.  6-4.  The  phase-only  reconstructions  after  10  iterations  using  (a)  Zf  =  7.3, 
(b)ZF  =  7.7,  and  (c)Zf  =  8  1 
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Figure  6-5  is  a  plot  of  the  errors  in  the  transform  plane,  Ep,  for  the  10 


Fig.  6-5.  Absolute  error  in  the  transform  plane,  Ep('>  vs.  the  number  of  iterations  using. 
The  assumed  number  of  Fresnel  zones  is  indicated. 

iterations  in  each  reconstruction.  Having  a  5%  error  in  the  number  of 
Fresnel  zones  is  seen  to  have  a  significant  effect  on  the  absolute  phase 
error  in  the  transform  plane.  The  two  iterations  using  the  incorrect  value 
for  the  number  of  Fresnel  zones  show  no  convergence  in  the  transform 
error. 

Figure  6-6  is  a  plot  of  the  relative  squared  errors  in  the  image  plane,  Ei, 
for  the  10  iterations  in  each  reconstruction.  Again,  using  the  incorrect 


2.5 


] 

I 

5 

i 


Zf  ■  8.1 


0.5 


7.7 

i  t  I  I _ I _ I _ 

4  5  6  7  8  8  10 

Itarotion  Nu«b«r 


Fig.  6-6  Relative  squared  errors  in  the  image  plane,  Ejl')  vs.  the  number  of  iterations 
The  assumed  number  of  Fresnel  zones  is  indicated. 

number  of  Fresnel  zones  has  a  significant  effect  on  the  error  in  the  image 
plane.  There  is  no  convergence  in  the  image  plane  error  as  the 
reconstruction  algorithm  progresses. 

6.2.3  Magnitude-only  reconstructions 

To  examine  the  sensitivity  of  the  magnitude-only  reconstruction 
algorithm,  a  set  reconstructions  using  magnitude-only  data  from  7.7 
Fresnel  zones  were  done  with  the  parameter  Zp  set  equal  to  7.3,  7.7,  and 
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8.1  in  the  reconstruction  algorithm.  This  is  the  same  type  of 
reconstruction  that  was  discussed  in  the  previous  example  of  phase-only 
reconstructions.  Figures  6-7(a),  (b).  and  (c)  are  the  results  of  this 
magnitude-only  reconstruction  after  20  iterations.  There  is  still  a 
recognizable  image  reconstructed  but  the  quality  of  the  reconstruction  is 
noticably  degraded  from  that  which  is  obtained  when  the  correct  number 
of  Fresnel  zones  is  used. 

Figure  6-8  is  a  plot  of  the  absolute  magnitude  error  in  the  transform 
plane,  Em.  for  the  20  iterations  in  each  reconstruction.  In  this  case,  there 
appears  to  be  a  convergence  of  the  reconstruction  algorithm,  but  the 
errors  are  consistently  higher  when  the  incorrect  number  of  Fresnel  zones 
is  used. 

Figure  6-9  is  a  plot  of  the  relative  squared  error  in  the  image  plane,  Ei, 
for  the  20  iterations  in  each  reconstruction.  Again,  there  is  a  convergence 
in  the  error  when  there  is  an  error  in  the  number  of  Fresnel  zones,  but  the 
stagnation  level  of  the  error  is  consistently  higher  when  the  incorrect 
number  of  Fresnel  zones  is  used. 

From  these  reconstruction  experiments  it  has  been  shown  that  the 
reconstruction  algorithm  is  sensitive  to  a  5%  error  in  the  value  for  the 
number  of  Fresnel  zones,  but  that  the  reconstruction  algorithm  is  not 
unstable  with  this  type  of  error.  In  reconstructing  from  phase-only  data, 
there  is  a  marked  degradation  in  the  reconstructed  image  as  the 
algorithm  progresses.  The  reconstructed  image  during  the  first  few 
iterations  is  fairly  insensitive  to  an  error  in  the  number  of  Fresnel  zones, 
but  as  the  reconstruction  algorithm  continues,  the  quality  of  the 
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reconstructed  image  is  sacrificed  in  an  attempt  to  minimize  the  overall 
error  in  both  the  image  and  transform  planes.  In  reconstructing  from 
magnitude-only  data,  there  is  a  definite  degradation  in  the  quality  of  the 
reconstructed  images  and  the  stagnation  levels  of  the  errors  are  larger 
when  there  is  an  error  in  the  number  of  Fresnel  zones.  The  effects  for  very 
large  errors  are  discussed  below  in  Sect.  6-3. 


(a)ZF  *  7.3 


(b)ZF  =  7.7 


(c)Zf  =  8.1 
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Fig.  6-7.  The  magnitude-only  reconstructions  after  20  iterations  using  (a)  Zf  =  7.3, 
(b)Zc  =  7.7,  and  (c)Zf  =  8  1. 
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Fig.  6-8.  Absolute  error  in  the  transform  plane.  £^1')  vs.  the  number  of  iterations .  The 
assumed  number  of  Fresnel  zones  is  indicated. 
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6.3  Searching  Through  the  Fresnel  Region 

As  discussed  above,  if  the  number  of  Fresnel  zones  is  known  only 
approximately,  the  reconstruction  algorithm  is  stable,  but  yields  a 
reconstructed  image  which  degrades  slowly  and  steadily  depending  upon 
how  close  to  the  correct  value  of  Zp  the  algorithm  is  using.  If  no 
information  about  the  number  of  Fresnel  zones  is  available,  a  guess  must 
be  used. 

6.3.2  Searching  using  Magnitude-only  data 

Figure  6-10  shows  the  magnitude  data  of  a  Fresnel  zones  transform 


Fig.  6- 1 0.  Magnitude  of  a  Fresnel  zone  transform  at  an  unknown 
location  in  the  Fresnel  zone  region. 

from  an  unknown  location  in  the  Fresnel  zone  region.  It  is  assumed  that 
because  of  sampling  considerations,  the  number  of  Fresnel  zones  is  less 
than  or  equal  to  12.0. 
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The  first  stage  of  the  search  for  the  correct  number  of  Fresnel  zones 
involves  trying  to  reconstruct  from  various  locations  in  the  Fresnel  zone 
region.  For  a  set  of  values  0.5  12.0  with  an  increment  of  0.5,  the 

magnitude-only  reconstruction  algorithm  was  run  through  two  complete 
cycles  for  each  value  of  If.  Both  the  absolute  magnitude  error,  Em.  and 
the  relative  squared  error,  Ei,  were  tabulated  for  the  attempted 
reconstructions.  Figure  6-11  is  a  plot  of  the  absolute  error  in  the 


Fig.  6-11.  Plot  of  the  absolute  magnitude  error,  EmO),  vs.  guesses  for  the  number  of 
Fresnel  zones. 

transform  plane  after  the  second  iteration  verses  the  assumed  number  of 
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Fresnel  zones.  Note  that  this  curve  has  a  distinct  minimum  in  the  area 
around  8.0  to  8.5  Fresnel  zones. 

To  obtain  a  more  accurate  estimate  for  the  unknown  value  of  Zf,  a  fine 
scale  search  is  now  done  for  8.0^  Zp^  8.5  with  step  a  size  of  0.1.  The 
results  of  this  fine  scale  search  are  shown  in  Fig.  6-12.  Again,  note  that 


Fig.  6-12.  Plot  of  the  absolute  magnitude  error,  vs.  fine  scale  guesses  for  the 
number  of  Fresnel  zones. 

there  is  a  distinct  minimum  in  the  region  near  8.3  Fresnel  zones. 

Using  the  value  of  8.3  for  Zp,  a  number  of  iterations  were  now 
performed.  The  reconstructed  image  after  10  iterations  is  shown  in  Fig.  6- 
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13.  This  reconstructed  image  is  easily  recognized.  It  is  interesting  to  note 


Fig.  6-13.  Reconstruction  after  10  iterations  using  Zf  =  8.3.  The  true 
number  of  Fresnel  zones  is  8.32. 

that  the  unknown  magnitude  data  shown  in  Fig.  6-10  was  actually 
calculated  at  8.32  Fresnel  zones.  This  searching  technique  could  be 
continued  to  yield  a  more  accurate  estimate  for  the  number  of  Fresnel 
zones,  but  it  was  demonstrated  in  the  previous  section  that  the  quality  of 
the  reconstruction  is  dependent  upon  the  estimate  for  the  number  of 
Fresnel  zones.  The  time  and  effort  spent  on  this  searching  would  depend 
upon  the  quality  of  the  reconstruction  which  is  desired. 

6.3.3  Searching  Using  Phase-only  data 

Figure  6-14  shows  the  phase  data  of  a  Fresnel  zone  transform  from  an 
unknown  location  in  the  Fresnel  zone  region.  The  same  type  of  searching 
algorithm  that  was  discussed  in  the  previous  section  was  used  with  the 


Fig.  6-14.  Phase  of  a  Fresnel  2one  transform  at  an  unknown  location  in 
the  Fresnel  zone  region. 

same  assumption  that  the  number  of  Fresnel  zones  is  less  than  or  equal  to 
12.0. 

The  first  stage  of  the  search  for  the  correct  number  of  Fresnel  zones 
requires  a  finer  search  through  the  entire  possible  range  of  Fresnel  zones. 
For  a  set  of  values  0.25  =Zf=  12.0  with  an  increment  of  0.25,  the  phase- 
only  reconstruction  algorithm  was  run  through  two  complete  cycles.  For 
each  guess  at  the  number  of  Fresnel  zones,  both  the  absolute  phase  error, 
Ep,  and  the  relative  squared  error,  Ei,  were  tabulated  for  the  for  the  first 
two  cycles  of  the  reconstruction  algorithm.  Figure  6-15  is  a  plot  of  the 
absolute  phase  error  in  the  transform  plane  after  the  second  iteration 
verses  the  assumed  number  of  Fresnel  zones.  Note  that  this  curve  has  a 
much  sharper  minimum  than  the  curve  shown  in  Fig.  6-11,  which  is  the 
reason  a  smaller  increment  in  the  number  of  Fresnel  zones  is  used  in  the 
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Fig.  6-15.  Plot  of  the  absolute  phase  error,  EpW,  vs.  guesses  for  the  number  of  Fresnel 
zones. 

searching.  The  very  sharp  minimum  is  in  the  area  around  8.0  to  8.5 
Fresnel  zones. 

To  obtain  a  more  accurate  estimate  for  the  unknown  value  of  Zp,  a  fine 
scale  search  is  now  done  for  7.8^  Zp^  8.7  with  step  a  size  of  0.1.  The 
results  of  this  fine  scale  search  are  shown  in  Fig.  6-16.  Again,  note  that 
there  is  a  distinct  minimum  when  the  guess  for  the  Fresnel  zones  is  equal 
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Fig.  6- 1 6.  Plot  of  the  absolute  phase  error,  Ep('),  vs.  fine  Kale  guesses  for  the  number  of 
Fresnel  zones. 

Using  the  value  of  8.3  for  Zp,  a  number  of  iterations  were  now 
performed.  The  reconstrurted  image  after  10  iterations  is  shown  in  Fig.  6- 
18.  This  reconstructed  image  is  easily  recognized.  It  is  interesting  to  note 
that  as  in  the  previous  example,  the  given  phase  data  shown  in  Fig.  6*14 
was  actually  calculated  at  8.32  Fresnel  zones.  An  exact  value  for  the 
correct  number  of  Fresnel  zones  is  not  needed  to  successfully  recover  a 
good  quality  image. 
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Fig.  6-17.  Reconstruction  after  5  iterations  using  Zf  =  8.3.  The  true 
number  of  Fresnel  zones  is  8.32. 

6.3.4  Conclusions 

In  Sect  6.2  it  was  shown  that  a  5%  error  in  the  number  of  Fresnel  zones 
produces  a  definite  degradation  in  the  quality  of  the  reconstructions. 
Both  the  magnitude-only  and  phase-only  reconstruction  algorithms  are 
sensitive  to  this  type  of  error. 

The  stability  and  sensitivity  of  the  reconstruction  algorithms  allowed 
for  the  development  of  another  algorithm  which  is  designed  to  search 
through  the  Fresnel  zone  region  to  determine  the  correct  value  of  Zf  The 
error  metrics  used  to  do  the  searching  are  the  same  measures  of  error  that 
are  used  to  monitor  the  progression  of  the  reconstruction  algorithms.  It  is 
very  likely  possible  to  develop  improved  error  metrics  and  searching 
techniques  for  the  reconstruction  of  images  from  partial  data  at  an 
unknown  location  in  the  Fresnel  zone  region. 
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Appendix  A 


Aberrated  Wavefront  Illumination 

An  alternative  view  of  the  optical  Fresnel  zone  transform  described  in 
Sect.  2.2  will  be  presented  in  this  appendix.  The  image  reconstruction 
algorithms  discussed  in  Chapt.  2  assumed  that  the  object ,  which  is  real¬ 
valued  and  non-negative,  is  illuminated  by  a  coherent  plane  wave. 
Consider  the  optical  transform  system  shown  in  Fig.  A-1.  In  this  case,  the 


Fig.  A-1  Coherent  optical  Fourier  transform  system. 

object  is  known  to  be  real-valued  and  non-negative,  the  Fourier  transform 
plane  is  located  in  the  back  focal  plane  of  the  lens,  and  the  illuminating 
wavefront  has  an  amplitude  and  phase  described  by  the  function, 
a(x,y)exp{iW(x,y)}. 
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The  scalar  components  of  the  electric  field  in  the  front  and  back  focal 
planes  of  the  lens  are  related  through  a  Fourier  transform  pair.  If  the 
amplitude  transmittance  of  the  image  is  denoted  by  g(x,y).  and  the  scalar 
component  of  the  field  in  the  transform  plane  is  denoted  by  G(fx,  fy),  then 
the  Fourier  transform  relation  which  relates  the  fields  in  the  front  and 
back  focal  planes  of  the  lens  is  given  by 


f  f w.  »  -iinif  x*fy) 

J  j  ^Uj')aUjr)  ‘  *  dxdy, 

and  the  inverse  transform  is  given  by 


(A-1) 


1  »  f  f  +i2n(r  x-¥f  y) 

tfU.y)  =  — I  Gif  f)t  '  ^  dfdf 

aOt.y)  *  y  *  y 

where  the  spatial  frequency  variables  fx  and  fy  have  been  defined  as 


(A-2) 


(A-3) 


f  = —u/kF  and  f  =—vlkF. 

Notice  that  if  a{x,y) » 1 .0  and  W(x,y)  =  -nai(x2  ♦  y2),  then  Eqs.  (A-1)  and  (A- 
2)  reduce  to  the  Fresnel  zone  transform  pair  described  in  Chapt.  2. 

It  is  now  possible  to  apply  the  reconstruction  techniques  described  in 
Chapt.  2  to  the  case  of  a  generalized  illuminating  wavefront.  Consider  a 
phase  function  for  the  illuminating  wavefront  which  contains  the  5 
primary  aberrations;  the  phase  of  this  wavefront  can  be  written  as 


W(x.>)  =  2n{AU*  +  >*)*+  ByU*  +  / )  +  C(x“  +  3y’ ) 


+  Dix^  +  y^)  +Ey  +Fx)  , 


(A-5) 


where  A  is  the  number  of  waves  of  spherical.  B  is  the  number  of  waves  of 
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coma.  C  is  the  number  of  waves  of  astigmatism.  D  is  the  number  of  waves 
of  defOvUS.  E  is  the  number  of  waves  of  tilt  about  the  x-axis.  and  F  is  the 
number  of  waves  of  tilt  about  the  y-axis. 

The  original  256  X  256  X  8-bit  digitized  image  used  in  the  following 
reconstruction  experiments  is  shown  in  Fig.  A-2. 


wavefront  is  very  similar  to  the  information  one  would  obtain  from  a 
Fourier  transform  obtained  with  plane  wave  illumination. 

The  magnitude-only  and  phase-only  reconstructions  are  shown  in  Fig. 
A-4.  The  magnitude-only  reconstruction  after  10  iterations,  shown  in  Fig. 
A-4(a),  is  easily  recognized  and  is  about  the  same  quality  as  a  magnitude- 
only  reconstruction  after  5  iterations  from  Fresnel  zone  data  at  3.0  Fresnel 
zones.  The  phase-only  reconstruction  after  5  iterations,  shown  in  Fig.  A- 
4(b).  is  very  good.  In  the  phase-only  reconstruction  the  edges  are 
recovered  very  quickly,  and  as  the  iterations  continue  the  low  frequency 
areas  of  the  picture  are  slowly  improved. 

A  second  set  of  reconstruction  experiments  was  performed  with  the 
parameters  a(x,y)=1.0,  A  =  4.0,  B=3.0,  C  =  -4.0,  D  =  -1.0,  E  =  0.0,  and 
F  =  0.0.  The  magnitude  and  phase  of  of  this  Fourier  transform  are  shown 
in  Fig  A-5.  Notice  that  the  magrtitude  information  shown  in  Fig.  A-5(a)  Is 
more  spread  out  than  in  the  case  of  a  Fourier  transform  obtained  with 
plane  wave  illumination,  and  that  there  is  no  longer  any  symmetry  to  the 
magnitude  data.  The  phase  information  shown  in  Fig.  A-5(b)  has  small 
regions  which  show  a  slowly  varying  structure  in  the  phase  data. 

The  two  reconstruction  from  this  data  are  shown  in  Fig  A-6.  After  10 
iterations  the  magnitude-only  reconstruction  shown  in  Fig.  A-6(a)  is  fairly 
good,  and  is  about  the  same  quality  as  a  magnitude-only 
reconstructionafter  5  iterations  from  Fresnel  zone  data  at  5.0  Fresnel 
zones.  The  phase-only  reconstruction  after  5  iteration  shown  in  Fig.  A- 
6(b)  is  very  good. 
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(a)  Magnitude 


(b)  Phase 


Fig.A-3  (a)magnitude  and  (b)  phase  of  Fourier  transform  using  an  aberrated 
wavefront  with  a(x,  y)  =  0.1,  A  =  0.4,  B  =  0.3,  C  =  -0.4,  D  =  -0.1,  E  =  0.0,  and 
F  =  0.0. 
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(b)  Phase-only  Reconstruction 


Fig.  A-4  (a)Magnitude-only  reconstruction  after  ten  iterations.  (b)Phase-only 
reconstruction  after  five  iterations. 
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(a)  Magnitude 


(b)  Phase 

(a)magnitude  and  (b)  phase  of  Fourier  transform  using  an  aberrated 
wavefront  with  a(x,  y)=  1.0,  A  =  4  0,  B  =  3.0,  C  =  -4.0,  D  =  -1.0,  E  =  0.0,  and 
F  =  0.0. 
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A  coherent  optical  system  which  has  an  aberrated  wavefront 
illuminating  a  real-valued,  non-negative  object  has  been  discussed  in  this 
appendix.  Two  examples  of  the  resulting  magnitude  and  phase  of  the 
Fourier  transform  have  been  presented.  Through  these  examples,  it  was 
demonstrated  that  it  is  possible  to  recover  an  image  in  just  a  few 
iterations  from  either  the  phase-only  or  magnitude-only  part  of  the 
Fourier  transform,  if  the  form  of  the  aberrated  illuminating  wavefront  is 
known.  For  a  set  number  of  iterations  there  is  an  increase  in  the  quality  of 
the  reconstructed  image  obtained  from  magnitude-only  data  as  the 
aberrations  on  the  illuminating  wavefront  become  more  severe.  In 
reconstructing  from  phase-only  data,  a  good  reconstruction  is  obtained 
which  is  fairly  insensitive  to  the  amount  of  aberration  on  the  illuminating 
wavefront. 

In  reconstructing  an  image  from  Fresnel  zone  data,  it  is  assumed  that 
the  unknown  image  is  real-valued,  non-negative,  has  a  known  size,  and 
that  the  given  data  comes  from  a  known  region  in  the  Fresnel  zone  of  a 
coherent  optical  processor.  In  reconstructing  an  image  from  partial 
Fourier  transform  data  obtained  from  a  coherent  optical  system  with  a 
known  aberrated  illumination  wavefront,  the  unknown  image  is  again 
assumed  to  be  real-valued,  non-negative,  and  to  have  a  known  size.  This 
unknown  image  corresponds  to  the  magnitude  of  the  wavefront  in  the 
front  focal  plane  of  a  coherent  optical  processor.  The  Fresnel  zone 
transform  discussed  throughout  the  body  of  this  dissertation  is  a  special 
case  of  the  Fourier  transform  relationship,  obtained  with  an  aberrated 
wavefront,  discussed  in  this  appendix. 
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